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ABSTRACT 

The  objective  of  this  thesis  is  to  show  that  the 
exact  solution  to  the  one-dimensional  Schrodinger  equa¬ 
tion  may  be  obtained  everywhere  in  terms  of  "smoothly- 
varying"  functions  at  a  speed  controlled  by  Hamilton's 
equations  if  the  numerical  techniques  employed  embody 
a  methodological  compatibility  with  the  conceptual  link 
between  wave  and  particle  mechanics. 

The  thesis  begins  by  examining  the  historical 
development  of  numerical  techniques  for  solving  the 
one-dimensional  Schrodinger  equation  in  the  context 
of  preserving  a  smooth  numerical  link  between  classical 
and  quantum  mechanical  solutions  of  equations  of  motion. 
The  problem  of  calculating  the  wave  function  in  the 
allowed  and  forbidden  regions  is  reformulated  as  that 
of  obtaining  a  correction  to  a  semiclassical  asymptotic 
form;  at  the  same  time,  this  introduces  the  need  for 
solving  a  Riccati  equation  for  a  function  <j)  (x)  and  an 
associated  "phase  correction  function"  f (x) .  The 
relation  between  these  functions  is  discussed  stressing 
the  properties  of  bounding,  the  behaviour  at  the  turn¬ 
ing  point  and  the  particular  solution  to  the  equation 
for  (p  (x)  which  satisfies  the  requirement  that  the  wave 
function  become  the  complex  exponential  of  a  modified 
form  of  Hamilton  characteristic  function.  A  new  method 
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for  obtaining  the  solution  to  a  particular  class  of 
Riccati  equations  (including  the  desired  solution  to 
$  (x) )  is  presented  in  the  context  of  a  "linearized 
iterative  sequence."  Its  implementation  in  two  rapid 
numerical  algorithms  is  discussed  and  tested.  Finally, 
results  are  given  which  conclusively  show  that  the  wave 
function  may  be  precisely  obtained  with  the  speed  of 
solving  the  Hamilton- Jacobi  equations  of  motion  in  the 
region  where  JWKB  approximants  are  valid. 
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CHAPTER  I 


INTRODUCTION 


An  important  mathematical  problem  in  quantum 
mechanics  is  that  of  solving  the  one-dimensional 
Schrodinger  equation 


+  k2 (x)  y (x)  =  0 


(1.1) 


where 


k2(x)  =  %  [E  -  V  (x)  ]  . 

fi 


(1.2) 


This  equation,  a  form  of  Sturm-Liouville  problem,  was 
first  investigated  in  a  purely  mathematical  context 
by  Carlini  (1817)  ,  Liouville  (1837)  ,  and  Green  (1837)  . 
The  advent  of  quantum  mechanics  revived  interest  in 
the  problem  and  the  search  for  useful  approximate 
solutions  led  to  the  rediscovery  or  development  of 
the  JWKB  approximation  (after  Jeffreys ^  (1923)), 

Wentzel,  Kramers,  and  Brillouin  (1926)), 


x 


y  (x) 


(  k(x)  ) 


exp[±i  (  k(x')dx'  +  e  (x) ]  . 


(1.3) 


This  approximation  is  valid  (that  is,  the  correction 
function  e (x)  is  negligibly  small)  only  in  an  asymp¬ 
totic  domain  where  k (x)  is  non-zero  and  [k1 (x)/k(x)] 
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is  sufficiently  small.  Zeros  of  k (x)  are  called 
"turning  points/'  a  name  which  alludes  to  the  pro¬ 
perties  of  the  corresponding  classical  motion.  The 
solution  to  Schrodinger 1 s  equation  in  the  form  (1.3) 
is  characterized  by  a  "phase"  describing  the  wave 
properties  of  the  probability  amplitude  for  a  quantum 
mechanical  particle  and  an  amplitude  describing  an 
essentially  classical  variation  of  probability  density 
with  the  reciprocal  of  the  particle  momentum.  The 
JWKB  approximation  is  of  interest  not  only  physically , 
because  of  the  interpretation  which  can  be  given  to  the 
amplitude  and  the  leading  term  of  the  phase,  but 
mathematically ,  since  these  quantities,  unlike  the 
wave  function  y(x),  do  not  oscillate  but  vary  smoothly 
in  a  fashion  controlled  by  the  classical  momentum. 

When  sufficiently  accurate,  the  JWKB  asymptotic  form 
is  thus  computationally  far  superior  to  direct  inte¬ 
gration  of  the  Schrodinger  equation  (1.1),  particularly 
when  the  de  Broglie  wavelength  2ir/k(x)  is  small. 

These  properties  of  the  JWKB  approximation  have 
stimulated  other  workers  to  investigate  its  properties 
further  with  the  intent  of  extending  its  range  of 
validity  and  accuracy  and,  possibly,  its  physical 
interpretation.  While  certain  aspects  of  the  mathe¬ 
matical  problem  of  such  asymptotic  series  solutions  to 
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these  equations  go  back  as  far  as  the  work  of  Green 

(1837),  discussions  having  a  direct  bearing  on  the 

physical  interpretation  of  solutions  to  Schrodinger 1 s 

equation  are  more  recent.  In  this  connection,  we 

shall  mention  the  efforts  of  two  workers  as  a  resume  of 

relevant  work. 

(2 ) 

Langer  considered  the  Schrodinger  wave  func¬ 

tion  to  have  the  form 

x 


y(x) 


exp  [  A 


a  (x ' 


A)  dx1] 


(1.4) 


from  which  one  obtains  a  Riccati  equation  for  o  (where 
A  is  large) 


g  +  a2  (x)  +  p2  (x )  =  0 


where 


(1.5) 


p(x)  =  [k(x)/A]  .  (1.6) 

If  a(x,A)  is  represented  as  a  series  expansion  in  A 
that  is 


00  a  .  (x) 

a(x,A)  =  l  [-^ — ] 

j=0  A J 


(1.7) 


the  differential  equation  becomes 
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(cJq  (x)  +  p2  (x)  )  + 


3 


+  l  o  Ax)  o._  (x)  ]  =0  . 
£=0  J  * 


(1.8) 


By  choosing  the  coefficients  a ^  (x)  so  that  individual 
terms  vanish,  that  is, 


ao(x)  =  ±ip  (x)  ,  aj  (x)  =  [-1/2  aQ  (x)  ]  (x) 


(1.9) 


the  equation  is  formally  solved.  Langer  saw  that  a 
solution  of  this  type  is  apparently  asymptotic  and  may 
not  be  extended  across  the  turning  points.  Since  a 


2 

zero  of  k  (x)  is  a  singular  point  of  the  asymptotic 


equation  exactly  solved  by  (1.3)  (with  e (x)  =  0) ,  a 

"Stokes'  phenomenon"  occurs,  that  is,  an  analytic 

solution  y  (x)  to  equation  (1.1)  which  is  asymptotically 

well  represented  by  a  particular  linear  combination  of 

the  asymptotic  forms  (1.3)  on  one  side  of  a  zero  of 
2 

k  (x)  will  not  be  asymptotically  well  represented  by 
the  same  linear  combination  of  the  corresponding  analy¬ 
tic  continuations  of  the  forms  (1.3)  on  the  other  side 
of  the  zero;  a  "Stokes'  jump"  occurs  in  the  linear 
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combination  coefficients.  Langer  showed  that  a 

"connection  formula"  fixing  the  Stokes'  jumps  can 

be  obtained  using  the  asymptotic  expansions  of  solu- 

2 

tions  to  Airy's  equation  (assuming  that  k  (x)  is 
approximately  linear  in  x  in  the  neighbourhood  of 
its  zeros) . 

There  are  several  known  alternative  expansions 

of  the  JWKB  forms  to  that  of  Langer.  They  all  agree 

up  to  second  order,  which  is  the  form  (1.3)  with  e  =  0. 

This  fact,  plus  the  difficulty  of  computing  higher 

terms,  their  lack  of  obvious  physical  interpretation, 

and  the  relatively  poor  subsequent  convergence  of  all 

such  expansions  has  led  to  the  widespread  adoption  of 

(1.3)  as  the  standard  JWKB  form. 

( 3 ) 

Kemble  ,  following  Zwaan,  considered  the  wave 
function  to  be  a  linear  combination  of  such  second- 
order  JWKB  forms,  that  is 

y(x)  =  a+W+(x)  +  a_W_ (x)  (1.10) 

in  the  classically  allowed  region  (where  k(x)  is  real) 
and 

y(x)  =  $+W+(x)  +  3_W_  (x)  (1.11) 

in  the  classically  forbidden  region  (where  k(x)  is 

/s 

imaginary)  with  W+ (x)  defined  by 
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W+(x) 


(  k(x) 


exp [ ±i 


x 

» 

k (x ' ) dx ' ]  . 


(1.12) 


Owing  to  the  Stokes  phenomenon,  a+ ,  a_  are  not  equal 
to  3  ,  3_;  by  developing  exact  equations  for  the  way 
in  which  these  asymptotically  constant  coefficients 
change  along  a  path  in  the  complex  plane  circum¬ 
navigating  the  turning  point  on  the  real  axis ,  Kemble 
succeeded  in  obtaining  a  rigorous  derivation  of  the 
connection  formulae  without  appealing  to  Airy's  equa¬ 
tions  as  a  model.  Kemble's  approach  is  of  interest  to 
our  work  because  of  its  emphasis  on  the  connection  to 
exact  solutions  to  (1.1)  of  the  coefficients  a+ (x) . 

The  work  of  Langer  and  Kemble  is  a  pinnacle  of 

achievement  in  the  formulation  of  the  JWKB  method 

although  the  literature  on  the  subject  before  and 

after  their  work  is  enormous.  A  much  more  recent  book 

( 4 ) 

by  Froman  and  Froman  (1965)  compiles  and  summarizes 
the  known  literature  and  the  reader  is  referred  to  this 
work  for  a  more  complete  bibliography  than  can  be  given 
here . 

One  further  analytic  outgrowth  of  the  JWKB 
method  is  of  direct  interest  to  us.  Regular  solutions 
of  Schrodinger ' s  equation  associated  with  scattering 
by  a  central  potential  are  asymptotically  characterized 
by  a  JWKB  solution  (that  is,  a  linear  combination  of 
forms  (1.3))  with  a  characteristic  amplitude  and 
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"phase-shift,"  the  latter  constant  fixing  the  ratio  of 

the  linear  combination  coefficients.  It  then  seems 

natural  to  characterize  the  exact  regular  solution 

everywhere  in  the  classically  allowed  region  by  a  JWKB 

form  connected  by  variable  phase  and  amplitude  functions 

which  would  assume  the  correct  asymptotic  values  at 

infinity.  Much  of  the  work  on  these  techniques  is  due 

(5 ) 

to  or  has  been  summarized  in  the  book  by  Calogero  . 
However,  all  of  Calogero' s  methods  have  a  serious 
defect:  the  differential  equation  for  the  variable 

phase  shift  and  all  subsequent  modifications  of  it 
possess  a  "step-like"  behaviour  in  the  phase  function 
(associated  with  trigonometric  cycles  in  the  wave  func¬ 
tion)  so  that,  although  this  correction  to  the  JWKB 
form  is  small,  its  rapid  variation  at  such  "steps" 
makes  the  technique  of  relatively  little  numerical 
advantage . 

The  advent  of  the  high-speed  digital  computer 
has  fostered  the  growth  of  numerical  techniques  in 
solving  the  one-dimensional  Schrodinger  equation. 
Essentially,  these  are  truncated  Taylor  series  expan¬ 
sions  which  are  used  recursively  to  calculate  the  wave 
function  and  its  derivatives  at  a  set  of  points  which 
are  separated  by  some  small  but  finite  length  or  step- 
Such  techniques  have  two  serious  limications. 
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One  is  that  truncated  Taylor  expansions  contain  an 
error  term  which  is  bounded  by  the  stepsize.  There¬ 
fore,  to  make  the  calculation  reliable,  it  is 
necessary  to  select  the  stepsize  to  be  significantly 
smaller  than  the  local  wavelength  2-n-/k(x)  of  the 
Schrodinger  equation.  The  second  limitation  is  due 
to  the  finite  precision  of  computers.  Since  two 
linearly  independent  solutions  to  Schrodinger ' s 
equation  exist,  imprecise  calculations  will  introduce 
a  component  of  the  irregular  solution  which  will  even¬ 
tually  diverge.  Thus,  in  making  the  stepsize  suffi¬ 
ciently  small  to  ensure  the  accuracy  of  the  truncated 
Tayler  series,  the  number  of  calculations  will  decrease 
the  precision  of  the  calculations  thereby  promoting  the 
growth  of  the  undesirable  and  non-physical  solution. 


In  an  effort  to  minimize  the  growth  of  the 

( 6 ) 

irregular  solution,  Gordon v  '  has  proposed  that  the 


potential  energy  be  regarded  as  a  collection  of  piece- 
wise  continuous  line  segments.  Then,  in  any  one  of 
these  segments,  Schrodinger ' s  equation  takes  the  form 


(1.13) 


dz 


which  has  Airy  function  solutions,  namely 


-1/3  v  3 


,2  3/2 

(t  z 


(1.14) 


< 


By  suitably  matching  these  wave  functions  at  the 
boundaries  of  the  line  segments,  the  wave  function 
for  the  entire  range  may  be  constructed.  Gordon's 
method  is  superior  to  a  direct  numerical  integration 
in  terms  of  the  accuracy  and  the  speed  with  which  the 
calculations  may  be  performed.  Unfortunately,  both 
techniques  are  numerical  methods  only  and  offer  no 
physical  insight  into  the  problem. 

From  the  preceding  survey  of  numerical  techniques 
for  solving  the  one-dimensional  Schrodinger  equation, 
one  obtains  a  feeling  of  methodological  discontinuity; 
in  the  presently  existing  repertory,  there  seems  to  be 
no  technique  which  takes  advantage  of  any  conceptual 
link  between  quantum  mechanics  and  classical  mechanics. 
In  classical  mechanics,  one  integrates  Hamilton's 
equations  or,  equivalently,  the  Hamilton-Jacobi  equa¬ 
tions,  and  the  stepsize  for  numerical  integration  is 
limited  only  by  the  rate  at  which  the  potential  func¬ 
tion  changes  with  respect  to  the  local  energy;  in  the 
exact  integration  of  the  Schrodinger  equation,  on  the 
other  hand,  the  de  Broglie  wavelength  2ir/k  (x)  determine 
the  stepsize  in  place  of  the  natural  stepsize  for 
Hamilton's  equations,  except  when  the  latter  is  smaller 
(as  in  the  case  of  low  energy) .  In  the  high  energy 
limit,  the  JWKB  approximation  recovers  an  estimate  of 


.  ' 
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the  solution  which  can  be  made  at  the  speed  of  Hamil¬ 
ton's  equations,  but  there  is  no  simple  technique 
which  permits  a  smooth  extension  by  numerical  methods 
to  the  exact  solution  without  introducing  the  dis¬ 
continuity  of  method  and  corresponding  difficulty 
represented,  for  example,  by  the  much  smaller  stepsize 
of  the  direct  integration  of  Schrodinger ' s  equation. 

The  motivation  for  this  work  lay  in  the  persistent 
feeling  that  this  discontinuity  of  numerical  method 
and  accompanying  numerical  difficulty  could  in  fact  be 
eliminated,  and,  as  we  shall  see,  this  feeling  is 
correct. 

A  clue  to  the  origins  of  the  discontinuity  of 
method  is  provided  if  one  examines  the  nature  of  the 
conceptual  link  between  quantum  and  classical  mechanics, 
and  this  link,  of  course,  involves  the  JWKB  approxima¬ 
tion.  Quantum  mechanics  replaces  the  classical  concept 
of  a  point  particle  in  motion  with  that  of  a  moving 
wave  packet.  The  size  of  the  packet,  a  measure  of  the 
uncertainty  in  the  particle's  position,  is  proportional 
to  Planck's  constant.  Thus,  in  the  limit  that  h  goes 
to  zero,  quantum  mechanics  and  classical  mechanics 

should  provide  the  same  description  of  matter.  It  is 

( 7 )  .  .  ... 

well-known  that,  in  this  limit,  the  Schrodinger 

equation  becomes  equivalent  to  the  Hamilton- Jacobi 
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equations  of  motion  provided  that  the  wave  function  has 
the  form 

y  (x)  =  exp  {iS(x)/h}  (1.15) 


where  S (x)  is  Hamilton's  characteristic  function, 
namely 


S(x) 


x 


/2m [E  -  V (x ' ) ]  dx ' 


(1.16) 


It  is  important  to  note  that  this  solution  is  identical 
to  Langer's  when  A  becomes  infinite.  (Langer's  defini¬ 
tion  of  <jq  (x)  includes  an  indeterminacy  in  its  sign 
as  a  result  of  taking  the  square  root  of  a  quantity. 

In  the  classical  picture,  this  corresponds  to  a  parti¬ 
cle  whose  momentum  has  a  well-defined  magnitude  but 
whose  direction  of  motion  is  purely  left  going  or 
purely  right  going.)  Moreover,  if  one  assumes  A  to 
be  1/fi,  then  Langer's  formal  solution  to  the  Riccati 
equation  may  be  interpreted  as  the  classical  solution 
followed  by  a  perturbation  series  in  fi.  Furthermore, 
the  term  p (x)  is  identically  the  classical  momentum. 
Therefore,  [-ia(x)]  is  the  "quantum  mechanical  momentum" 
and  possesses  the  classical  momentum  as  its  leading 
term  followed  by  a  perturbation  series  in  "h.  Indeed, 
this  is  consistent  with  the  momentum-operator  equation 


' 
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-ifi 


dy  (x) 
dx 


[-i o  ( x )  ]  y  (x) 


(1.17) 


If,  therefore,  one  wishes  to  preserve  a  smooth 
numerical  link  between  classical  and  quantum  mechani¬ 
cal  solutions  of  equations  of  motion,  the  preceding 
discussion  strongly  suggests  that  one  pay  special 
attention  not  to  the  general  solution  of  Schrodinger ' s 
equation  but  to  those  two  particular  solutions  which, 
in  the  classical  limit,  corresDond  to  a  momentum 
eigenfunction,  that  is  the  solutions  of  Schrodinger 1 s 
equations  which  also  asymptotically  obey  (1.17).  We 
shall  see  that  the  greater  numerical  difficulty  of 
previous  techniques  for  exact  solution  of  Schrodinger ' s 
equation  is  entirely  the  result  of  attempting  a  direct 
calculation  of  a  general  solution  of  the  form 

x 


y  (x) 


A  (x) 


sin  r 
cos 1 

✓ 


p  (x1  )  dx  '  +  6  (x)  } 


(1.18) 


Such  a  wave  function  poses  conceptual  difficulties  when 
one  tries  to  make  the  transition  to  classical  mechanics, 
because  it  describes  the  superposition  of  two  equi- 
probable  states  of  a  quasi-classical  particle  (left 
and  right  going  waves);  conversely,  the  apparently 
discontinuous  increase  in  the  difficulty  of  obtaining 
an  exact  numerical  solution  on  going  from  classical  to 


13 


quantum  mechanics  might  be  due  entirely  to  using  (1.18) 
rather  than  the  "complex  exponential"  solutions  as  the 
particular  exact  solutions  to  be  generated.  In  other 
words ,  the  intimate  asymptotic  connection  between 
complex  exponential  solutions  and  Hamilton's  charac¬ 
teristic  function  suggests  that  they  can  be  exactly 
calculated  with  the  ease  and  speed  of  solving  the 
Hamilton-Jacobi  equations. 
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CHAPTER  II 

REFORMULATION  OF  THE  PROBLEM  AS  A  CORRECTION  TO  AN 

ASYMPTOTIC  FORM 


A.  Classically  Allowed  Region 

Define  functions  a+  (x)  so  that  the  solution  to 
the  one-dimensional  Schrodinger  equation  has  the  form 


y(x)  =  a+(x)W+(x)  +  a_(x)W_(x) 

(2.1) 

y '  (x)  =  a ,  (x)  W  (x)  +  a  (x)W  (x) 


where 


W±  (x) 


x 

» 

exp  {±i 


k  (x  '  ) dx '  } 


(2.2) 


This  solution  corresponds  with  a  linear  combination  of 
two  asymptotic  exponential  solutions  to  Schrodinger ' s 
equation.  The  coefficients  are  allowed  to  vary  con¬ 
tinuously  so  that  the  wave  equation  will  be  satisfied 
exactly.  Our  selection  of  W+ (x)  to  be  non-singular  at 
the  turning  point  (unlike  Kemble's  choice)  suggests 
that,  by  suitably  obtaining  the  coefficients  a+ (x) , 
the  wave  function  will  be  regular  everywhere. 

By  simultaneously  solving  equations  (2.1),  the 


coefficients  are  observed  to  be 


. 
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a+  (x) 


I 


(x)W_ 

± 


(x)  -  y  (x)  (x) 

-2ik (x) 


(2.4) 


This  expression  may  be  simplified  by  defining  a  new 
set  of  coefficients  3  +  (x)  by 


3+(x)  =  a+(x)  {-2ik  (x) } 


(2.5) 


so  that 


3  ,  (x)  =  ±  (y(x)W_(x)  -  y  (x)W_(x)} 

+  + 


(2.6) 


Differentiating  once  and  employing  the  explicit  form 
for  W+ (x)  and  the  wave  equation  for  y(x),  the  problem 
reduces  to  that  of  solving 


3+ (x)  =  g(x)(3+(x)  +  3_(x)exp[  + 

+ 


h (x) dx]  } 


x 


(2.7) 


where  g  (x)  and  h (x)  are  defined  by 


g  (x)  = 


k'  (x) 
2k  (x) 


(2.8) 


h  (x)  =  2ik(x) 


This  equation  may  be  further  simplified  by  eliminating 
the  diagonal  coupling  in  the  equation.  Defining  new 
coefficients  y+  (x)  by 


16 


x 


3.  (x)  =  y ,  (x)  exp{ 


g  (x  ' ) dx '  } 


one  finds  that 


(2.9) 


Y_  (x)  g  (x) 
+ 


exp{  + 


x 


h (x ' ) dx ' }  . 


(2.10) 


Coupled  first-order  differential  equations  of  this  form 
have  been  investigated  in  the  theory  of  inelastic 
atomic  collisions.  In  particular,  Lebeda  and  Thor s on v  y 
found  that  undesirable  and  spurious  double-frequency 
oscillations  in  numerical  solutions  to  such  equations 
were  generated  by  inappropriate  usage  of  initial  condi¬ 
tions  and  that  these  could  be  suppressed  by  introducing 
a  suitable  asymptotic  solution  form.  They  defined  co¬ 
efficients  b+ (x)  and  b_ (x)  by 


b+(x)  =  y+(x) 


x 


(2.11) 


b_ (x)  =  y_ (x)  exp{- 


h (x 1 ) dx ' } 


from  which  one  obtains  the  following  pair  of  coupled 
equations : 


b_(x)  +  h(x)b_(x)  =  b+(x)g(x) 


b+(x)  =  b_(x)g(x) 


(2.12) 


' 
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Lebeda  and  Thorson  were  led  to  introduce  such  an  asym¬ 
metric  bias  into  the  symmetric  equations  (2.10)  by  the 
physical  initial  conditions  for  inelastic  scattering 
which  provided  that  asymptotically  y+(x),  say,  has 
the  value  1  while  y_ (x)  goes  to  zero  (in  that  problem, 
the  coefficients  represent  amplitudes  for  initial  and 
final  states,  respectively,  before  collision).  Why 
should  solutions  possessing  an  asymptotic  asymmetrical 
bias  of  this  sort  be  of  interest  in  our  present  problem? 
A  posteriori  we  may  appeal  to  the  intuition  of  Chapter  I 
which  emphasizes  that  it  is  the  "pure  complex  exponen¬ 
tial"  solutions  which  bear  a  direct  interpretive  connec¬ 
tion  to  the  motion  of  a  classical  particle. 

Following  Lebeda  and  Thorson,  therefore,  suppose 
that  we  consider  a  particular  solution  of  (2.1)  for 
which  a_ (x)  is  very  much  smaller,  at  least  asymptoti¬ 
cally,  than  a+ (x) ;  then  y_  (x)  will  be  much  smaller 
than  y+ (x) .  Since,  then,  y_ (x)  will  be  largely  "driven" 
by  y+ (x)  via  the  coupled  equations  (2.10),  it  is  reason¬ 
able  and  consistent  to  think  that  the  asymmetrically 
biased  quantities  b+ (x)  will  both  be  slowly  varying. 

The  coupled  equations  for  b+  (x)  may  be  converted  to  a 
single  Riccati  equation 


g(x) 

h  (x) 


1  (x) 
h  (x) 


(x) 


{1  -  <p2  (x)  } 


(2.13) 
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b+(x)  =b+(xQ)exp{  (j>  (x'  )g(x'  )dx'  }  , 


(2.14) 


where 


(x)  b+  (x)  =  b_  (x) 


(2.15) 


Now,  any  attempt  to  numerically  integrate  equation 
(2.13)  would  only  re-introduce  the  oscillations  we 
wish  to  avoid.  However,  it  is  immediately  of  interest 
to  observe  that  far  from  a  turning  point,  an  approxi¬ 
mate  solution  of  (2.13)  is 


(2.16) 


indeed,  if  [g(x)/h(x)]  is  sufficiently  small,  it  is 
evident  that  such  a  solution  for  cp  would  justify  the 
rationale  in  selecting  the  transformations  for  b+ (x) , 
namely  that  b_  (x)  is  small  and  b+ (x)  is  slowly  varying. 
But  [g(x)/h(x)]  small  is  just  the  condition  for  asymp¬ 
totic  validity  of  the  JWKB  approximation. 

It  is  now  desirable  to  construct  the  wave  func¬ 
tion  in  terms  of  (j)  (x)  .  It  is  useful  to  observe  that 


x 


exp{ 


(2.17) 


x 


o 


. 
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Combining  these  with  the  transformations  made  earlier, 
the  a+  (x)  coefficients  become 

x 


a+  (x) 


a+(xo) 


^  ^  j  exp{  <f>  (x1  )g  (x1  )dx'  } 

L  X 


(2.18) 


a_  (x)  =  (})  (x)  exp  { 


h (x ' ) dx' }a+ (x) 


Is  this  result  consistent  with  the  assumptions  made  in 
constructing  the  various  transformations,  namely  that 
the  magnitude  of  a_  (x)  is  much  smaller  than  the  mag¬ 
nitude  of  a+ (x) ?  Since  h(x)  is  pure  imaginary,  the 
exponential  term  does  not  alter  the  magnitude  of  a_  (x) . 
Thus,  only  the  magnitude  of  <{>  (x)  is  significant.  By 

equation  (2.15)  ,  (J)  (x)  can  be  approximated  by  k'  (x)/ 

2 

4ik  (x)  which  is  small  in  the  region  where  the  JWKB 
approximation  is  valid.  Therefore,  equations  (2.18) 
are  indeed  consistent  with  the  assumptions  that  have 
been  employed. 

With  this  understanding,  the  wave  function  and 
its  derivative  (neglecting  proportionality  constants 
common  to  both)  may  be  written 

x 


-V 


y  (x)  =  [1  +  <J>  (x)  ]  [k  (x)  ]  2  exp{ i 


[k  (x)  - 


ik 1 (x) 
2k  (x) 


<p  (x)  ]  dx} 

(2.19) 


y'(x)  =  i[l-  <j)  (x)  ]  [k  (x)  ]  exp{i 


[k  (x)  - 


ik' (x) 


X 


2k  (x) 


<J>  (x)  ]  dx}  . 


1 


’ 


‘ 


‘ 
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Note  that  when  <|>(x)  is  zero,  these  expressions  are 
identically  those  of  the  "pure  exponential"  second 
order  JWKB  approximant.  Thus,  even  the  crudest  estimate 
of  cf)  (x)  can  give  very  accurate  wave  functions.  This 
form  preserves  the  wave  character  of  the  solution 
(as  a  pure  right  going  wave)  and  describes  the  classi¬ 
cal  variation  of  probability  density  with  the  reciprocal 
of  the  particle's  momentum.  However,  the  nature  of  the 
solution  at  the  turning  point  is  unknown  and  the  general 
form  of  the  solution  is  not  conveniently  that  of  a  com¬ 
plex  exponential.  We  shall  discover  later  how  these 
difficulties  may  be  circumvented. 

In  the  preceding  discussion,  it  was  assumed  that 
a+(x)  dominated  a_(x),  that  is  to  say,  the  wave  was 
right  going.  Consider,  to  the  contrary,  that  a_ (x) 
dominates  a+(x);  i.e.  that  the  wave  is  left  going. 

Then,  select  coefficients  b+ (x)  and  b_  (x)  by 


b  (x)  =  y  (x)  exp  { 


h (x) dx  } 


x 


o 


b_  (x)  =  y_  (x) 


(2.20) 


with  the  corresponding  pair  of  coupled  equations 

b+(x)  -  h(x)b+(x)  =  b_(x)g(x) 

■ 

b_  (x)  =  b+  (x)  g  (x) 


(2.21) 


' 


1 


' 
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This  may  be  converted  to  a  single  Riccati  equation, 
namely 


(x) 


g(x) 

h(x) 


{1-  <j,2  (x)  }  + 


<i>'  (X) 

h  (x) 


(2.22) 


where 


(x) 


b+  (x) 
b  (~x) 


(2.23) 


An  approximate  solution  to  this  equation,  valid  far 
from  the  turning  point  is 


<j>  (x) 


g  (x) 
h  (x) 


(2.24) 


It  is  important  to  note  that  equation  (2.22)  is  the 
complex  conjugate  of  equation  (2.13);  therefore,  one 
need  only  calculate  (p  (x)  for  one  of  these  equations 
since  the  solution  to  the  other  will  be  its  complex 
conjugate.  For  this  reason,  cp  (x)  will  be  used  to  denote 
the  solution  to  (2.13)  whereas  <p*  (x)  will  denote  the 
solution  to  (2.24).  In  identical  fashion,  one  obtains 
the  wave  function  and  its  derivative  for  a  left  going 
wave  to  be 


y  (x) 


y'  Cx) 


[1  + 


ik*  (x) 
2k  (x) 


cf)  *(x)  ]  dx} 


(2.25) 


[1- 


ik1 (x) 
2k  (x) 


4>*  (x)  ]  dx}  . 


* 


' 


. 


■ 


. 
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This  solution  to  the  wave  function  may  immediately  be 
recognized  as  the  complex  conjugate  to  (2.19).  Thus, 
any  general  solution  to  Schrodinger ' s  equation  may  be 
expressed  as  a  linear  combination  of  (2.19)  and  its 
complex  conjugate;  it  is  significant  that  this  includes 
all  real  solutions.  This  technique  must  now  be  applied 
to  the  task  of  calculating  the  wave  function  in  the 
forbidden  region. 


B.  Classically  Forbidden  Region 

The  first  portion  of  this  chapter  formulated  a 
technique  for  obtaining  the  solution  to  the  one¬ 
dimensional  Schrodinger  equation  in  the  allowed  region 
as  a  correction  to  an  asymptotic  form.  To  extend  this 
technique  into  the  forbidden  region,  consider  the  solu¬ 
tion  to  the  wave  equation  to  have  the  form 


y(x)  =  a, (x)W, (x)  +  a  (x)W  (x) 

-r  t  —  — 

y '  (x)  =  a+  (x) W+  (x)  +  a_(x)W_(x) 


(2.26) 


where 


x 


W+  (x)  =  exp  { ± 


k (x) dx} 


(2.27) 


' 


* 
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The  detailed  calculations  are  identical  in  form  to 
those  employed  earlier.  For  the  sake  of  brevity,  one 
need  only  recognize  that  the  mapping 

k  (x)  - ►  -ik  (x)  (2.28) 

will  transform  the  equations  which  describe  the  allowed 
region  into  those  which  describe  the  forbidden  region. 
Thus,  g(x)  and  h(x)  become 


g  (x) 


k'  (x) 
2k  (x) 


h(x)  =  2k  (x) 


(2.29) 


For  the  case  that  a+  (x)  dominates  a_  (x) ,  the  solution 
to  the  wave  equation  may  be  written 


x 


y  (x)  =  [1+  <j>  (x)  ]  [k  (x)  ]  2exp{ 


[k(x)+^-~v  <j>(x)]dx} 


2k  (x) 


x 


(2.30) 


x 


y '  (x)  =  [1  -  <p  (x)  ]  [k  (x)  ]  ^  exp{ 


[k  (x)  +  cf>(x)]dx} 


X 


o 


where  <p  (x)  satisfies 


♦  W  =  U  -  *2(x)} 


<P 1  (x) 

h  (x) 


(2.31) 


Note  that,  in  the  region  far  from  the  turning  point, 
4>(x)  is  approximately 


■ 


. 


' 
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4>  (x) 


g  (x) 

h  (x) 


(2.32) 


For  the  case  that  a_ (x)  dominates  a+(x),  the  solution 
to  the  wave  equation  may  be  written 


— 


y  (x)  =  [1  +  <j>  (x)  ]  [k  (x)  ]  2  exp{ 


x 


+ 


y  '  (x)  =  [1  -  cj>  (x)  ]  [k  (x)  ]  2exp{ 


[-k  (x)  +  cf>  (x)  ]  dx} 

°  (2.33) 

[-k  (x)  +  <f>  (x)  ]  dx} 


x 


where  ({)  (x)  satisfies 


*(x)  =  -  a|4  u  -  $2(k)}  +  CJ4 

Y  h  (x)  '  h  (x) 


(2.34) 


Far  from  the  turning  point,  <p  (x)  is  given  approximately 
by 


<P  (x) 


g  (x) 

h  (x) 


(2.35) 


In  both  of  these  cases,  cp  (x)  is  a  real-valued 
function,  unlike  its  counterpart  in  the  allowed  region 
which  has  both  real  and  imaginary  components.  Thus, 
it  is  no  longer  possible  to  construct  a  solution  to 
the  latter  equation  for  (p  (x)  as  the  complex  conjugate 
of  the  solution  to  (2.31).  However,  this  is  not  impor¬ 
tant  since  only  one  solution  to  the  wave  equation  is 
required  (that  is,  the  decaying  mode)  in  the  forbidden 
region . 

■k  *  * 


' 
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This  chapter  commenced  by  reformulating  the  pro¬ 
blem  of  solving  the  one-dimensional  Schrodinger  equation 
as  a  correction  to  asymptotic  solutions  of  exponential 
form.  In  pursuing  this  approach,  a  pair  of  solutions 
was  obtained  in  both  allowed  and  forbidden  regions 
which  bear  a  striking  resemblance  to  pure  exponential 
JWKB  approximants ,  and  which  reduces  to  them  in  the 
asymptotic  limit  (k'(x)/k(x)  and  its  derivatives  suffi¬ 
ciently  small) .  Such  a  form  preserves  the  concept  of 
the  phase  of  a  wave  and  the  classical  variance  of  the 
probability  density  with  reciprocal  particle  momentum. 

In  the  semiclassical  domain,  a  good  approximate  solution 
for  the  function  <p  (x)  will  generate  a  solution  y(x) 
which  is  more  accurate  than  the  JWKB  solution. 

The  crucial  question,  of  course,  is  the  nature 
of  the  particular  solution  (p  (x)  to  equation  (2.13)  which 
is  embodied  in  the  solution  y  (x)  in  (2.19)  and  whose 
asymptotic  leading  term  is  the  small  quantity  given  by 
(2.16).  Is  this  cj)  (x)  bounded,  especially  at  the  turn¬ 
ing  point?  Is  it  a  smooth  function  of  x  or  does  it 
exhibit  oscillatory  or  step-like  behaviour  which  would 
destroy  any  numerical  utility  of  this  approach?  We 
were  originally  led  to  consider  the  analogy  with  the 
problem  of  Lebeda  and  Thorson  because  they  found  that 
an  iterative  solution  of  their  Riccati  equation  analogous 


’ 

1 

' 

. 

* 
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to  (2.13),  starting  with  the  zero-order  approximation 

(2.16),  was  indeed  a  smooth  function  without  the 

spurious  oscillations  they  wished  to  eliminate.  (The 

meaning  of  an  "iterative  solution"  to  a  non-linear 

equation  may  be  clarified  by  consulting  a  reference 

(9 ) 

such  as  Mathews  and  Walker  .  It  is  explicitly  shown 
there  that,  in  the  case  of  the  JWKB  problem,  iteration 
gives  rise  to  the  first  two  terms  of  the  Langer  series.) 
It  turns  out  that  a  modified  form  of  the  iterative 
method  generates  a  smooth  <f>  (x)  but  we  defer  discussion 
of  this  subject  until  the  iterative  solution  of  equa¬ 
tion  [(2.13),  (2.19),  (2.31),  (2.34)]  is  described  in 

Chapter  IV.  The  boundedness  of  cf>  (x)  is  significant 
since  we  hope  to  obtain  a  form  for  the  wave  function 
which  remains  non-singular  through  the  turning  point 
region.  Finally,  the  form  of  the  solution  (2.19) 
appears  to  abandon  our  original  suggestion  that  a 
"physical"  solution  should  be  the  pure  exponential 
of  a  modified  Hamilton's  characteristic  function,  that 
is,  a  smoothly-varying  function  of  the  potential.  In 
the  next  chapter,  we  shall  resolve  most  of  these 
questions . 


V 


J 
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CHAPTER  III 

PURE  EXPONENTIAL  WAVE  FUNCTION 


A.  Derivation  of  Pure  Exponential  Form  of  Solution 

The  preceding  chapter  succeeded  in  reformulating 
the  problem  of  solving  the  one-dimensional  Schrodinger 
equation  as  a  correction  to  an  asymptotic  form.  As  a 
consequence,  we  seemingly  have  abandoned  the  ideal  of 
obtaining  a  pure  exponential  solution.  Consider,  then, 
the  ratio  of  equations  (2.19), 


y'  (x) 
y  (x) 


ik  (x) 


rl  -  4>  (x)  i 
ll  +  <j>(x)' 


(3.1) 


This  may  be  regarded  as  a  linear  first  order  differen¬ 
tial  equation  in  y (x) .  The  corresponding  first  inte¬ 
gral  is 


y  (x) 


x 


exp  {i 


k  (x) 


fl  ~  4>(x) 

L1  +  <j>  (x) 


]  dx} 


(3.2) 


where  (p  (x)  is  determined  by  (2.13).  A  second  solution 
which  is  valid  in  the  allowed  region  is  its  complex 
conjugate 


x 

» 

k  (x) 


^  1  -  <P  *  ( x) 


]  dx}  . 


y* (x)  =  exp  {-i 


1  +<|>*(x) 


(3.3) 


' 
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Similarly,  the  solutions  which  are  valid  in  the  forbid¬ 
den  region  are 


y(x)  =  exp  {± 


k(x)  [ 


1  -  <J>±  (x) 
1  +  <p+  (x) 


•]  dx} 


(3.4) 


x 


where  cj>  (x)  and  (J)_  (x)  satisfy  (2.31)  and  (2.34)  res¬ 
pectively.  Now  that  solutions  to  the  wave  equation 
which  are  simultaneously  momentum  eigenfunctions  have 
been  formally  defined,  it  is  necessary  to  investigate 
the  properties  of  $  (x)  and  its  calculation. 


B.  Phase  Correction  Function 

The  pure  exponential  character  of  the  solution 
just  obtained  imparts  much  significance  to  a  quantity 
we  shall  call  f (x)  where 


f  (x) 


1  ~  4>  (x) 

l  +  4>  (x) 


(3.5) 


The  function  f  (x)  has  the  effect  of  modifying  the  JWKB 
phase;  for  this  reason,  it  shall  be  referred  to  as  the 
"phase  correction  function".  However,  unlike  the  phase 
correction  function  employed  by  some  techniques,  this 
one  provides  the  amplitude  of  the  wave  function  as  well. 

The  function  <p  (x)  is  related  to  f  (x)  by 

1  -  f(x) 

1  +  f  (x) 


<f>  (x) 


(3.6) 


' 


. 


■ 


. 
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Equations  (3.5)  and  (3.6)  are  examples  of  linear  frac¬ 
tional  transformations .  Such  a  transformation 
maps  lines  into  circles  and  circles  into  lines.  Of 
particular  interest  is  the  fact  that  (3.5)  maps  the 
region  of  complex  $- space  bounded  on  the  outside  by 
the  unit  circle  into  the  right-half  complex  plane. 

This  symmetrization  may  be  extended  to  the  dif¬ 
ferential  equations  which  define  <j>  (x)  ,  namely 


<f>  (x) 


G(x) 
2H  (x) 


{1  -  cf>2  (x)  } 


4)'  (x) 

H  (x) 


f  (x) 


H  (x) 
2G  (x) 


{1  -  f2 (x) } 


f  1  (x) 
G  (x) 


(3.7) 


where 


G(x) 


k'  (x) 
k(x) 


(3.8) 


and 


H  (x) 


2ik  (x) 
2k  (x) 


(allowed  region) 
(forbidden  region) 


(3.9) 


Therefore,  we  see  that  if  a  formal  solution  (p  (x)  be 
obtained  as  some  functional  Q [G (x)  ,H  (x) ] ,  then  the 
functional  Q[H(x),G(x)]  is  a  solution  for  f  (x) ; 
however,  it  is  not ,  in  general,  the  corresponding 
solution  defined  by  (3.6).  Hence,  such  a  symmetrical 


relation  is  not  of  direct  interest. 


■ 
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Certain  symmetry  properties  of  sets  of  solutions 
<J>  (x)  and  their  corresponding  f  (x)  are  of  interest.  We 
noted  earlier  that  the  complex  conjugate  of  any  solu¬ 
tion  to  equation  (2.13)  is  a  solution  to  equation 
(2.22).  In  addition,  the  reciprocal  of  any  solution 
to  (2.13)  (or  (2.31)  in  the  forbidden  region)  is  a 
solution  to  (2.22)  (or  (2.34)  in  the  forbidden  region). 
An  equally  intriguing  result  is  that,  if  <p  (x)  is  a 
solution  to  (2.31)  and  corresponds  with  some  f (x) , 
then  l/(p*  (x)  is  also  a  solution  to  (2.31)  which,  by 
inspection,  corresponds  with  -f*(x).  Extending  this 
unusual  relationship  between  the  equations  which  define 
4>  (x)  for  left  and  right  going  waves  to  the  phase  correc¬ 
tion  function,  this  corresponds  with  replacing  f  (x)  in, 
say,  the  left  going  solution 


y  (x)  = 


x 


exp  {i 


k(x)  f lef t (x) dx} 


exp  { 


k (x)  fleft(x) 


(allowed  region) 


(3.10) 


(forbidden  region) 


with  -f (x)  in  the  right  going  solution 

x 

exp  {-i 


y  (x)  = 


k(x)  f  r£ght  dx^  (aH°wed  region) 


x 


x 


exp  {- 


(3.11) 

k  (x)  f r j_g^t  (x)dx)  (forbidden  region) 


' 

■ 
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and,  consequently,  recovers  an  identical  solution  to 
the  wave  equation.  Yet  another  way  of  looking  at  this 
relationship  between  the  equations  for  <j>  (x)  is  that 
they  are  identical  up  to  the  selection  of  the  branch 
cut  in  evaluating  k  (x) .  Thus,  if  one  knows  the  general 
solution  to  cj)  (x)  determined  by  one  differential  equa¬ 
tion,  the  general  solution  to  the  complementary  dif¬ 
ferential  solution  is  also  known. 

The  initial  condition  for  the  equation  for  <j>  (x) 

in  the  allowed  region  was  observed,  in  the  preceding 

2 

chapter,  to  be  that  cf>  (x)  varies  as  k'  (x)/4ik  (x)  far 
from  the  turning  point.  Is  this  consistent  with  the 
pure  exponential  form  of  solution  obtained  in  this 
chapter?  Consider  the  integrand  of  (3.2) 

k(x)f  (x)  = 

=  k  (x)  {1  -  2<J>  (x)  +  0[4>2  (x)  ]  }  .  (3.12) 

Compare  this  with  the  first  two  terms  of  the  JWKB  (or 
Langer)  series 

k<x)  +  %rr!r  =  *<x)  •  (3.i3> 

ZK{X>  2kz  (x) 

2 

By  associating  -2<J>  (x)  with  +ik'(x)/2k  (x)  ,  the  com¬ 
parison  is  complete.  Moreover,  this  shows  that  f (x) 
is  given  approximately  by 


' 

'  i  , 
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f(x)  =  1  +  1-'-(x).  (3.14) 

2kz (x) 

far  from  the  turning  point.  Similar  results  hold  in 
the  forbidden  region. 

The  exponential  solution  which  has  been  obtained 
describes,  in  the  semiclassical  limit,  a  particle  mov- 
in  a  particular  direction.  The  phase  correction 
function  f  (x)  provides  for  deviation  in  amplitude  and 
phase  of  a  quantum  mechanical  particle  from  its  classi¬ 
cal  counterpart:  however,  if  the  qualitative  concept 
of  a  particular  direction  to  the  particle's  motion 
(and,  hence,  the  specific  idea  of  a  "pure  complex 
exponential"  solution)  is  to  retain  any  validity,  the 
complex  function  f (x)  must  remain  in  the  same  half  of 
the  complex  plane  (namely,  the  right  half)  for  all 
points  x  in  the  classically  allowed  domain.  From  the 
linear  fractional  transformation  (3.6),  this  in  turn 
requires  that  cj)  (x)  lie  within  the  unit  circle.  In 
Appendix  A,  a  proof  is  given  that  if  $  (xq)  satisfies 
this  condition  at  some  point  xq  in  the  classically 
allowed  domain,  then  $ (x)  will  be  so  bounded  for  all 
points  in  the  domain.  Therefore,  the  above  concepts 
remain  valid.  This  illustrates  a  crucial  difference 
between  the  problem  of  calculating  (x)  (or  f  (x)  )  and 
that  of  calculating  phase  and  amplitude  factors  asso¬ 
ciated  with  other  techniques:  the  function  <f>  (x)  is 


r  ‘d 


■ 


1 
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bounded!  All  that  is  now  required  is  a  technique  for 
calculating  it  in  various  regions. 

C .  Iterative  Solution  to  the  Riccati  Equation  for  $ (x) 

Consider  the  general  first  order  Ricatti  equa¬ 
tion 

i|;(x)  =  a  (x)  +  3(x)ip2(x)  +  y(x)i|;  (x)  .  (3.15) 

The  following  iterative  scheme  is  proposed  when¬ 
ever  cx(x),  3 (x)  and  y  (x)  are  analytic,  "small"  compared 
with  unity,  and  have  slowly  varying  derivatives. 

For  i  =  0,  1,  2,  ...,  a  solution  to  (3.15)  may 
be  obtained  by  solving  the  ancillary  equation 

tJk(x)  =  ai(x)  +  3i(x)ip?(x)  +  y^(x)t|k(x)  (3.16) 

by  invoking  a  natural  recurrence  scheme 

^(x)  =  +  ai(x)  (3.17) 

where 

6i  (x)  =  1  -  2ai  (x)  3i  (x) 

ai+l (x)  =  {3i(x)a?(x)  +  yi (x) ai (x) }/6i (x) 

(3.18) 

3i+i(x)  =  3i(x)/5i(x) 

Yi+1 (*)  =  Y±  (x)/6i (x) 


■ 
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given  that 


aQ  (x) 

=  a  (x) 

B0(x) 

=  3(x)  (3.19) 

Y0(x> 

=  y  (x) 

From  these  relations,  it  follows  that  the  required 
solution  to  (3.15)  is  given  by  (x)  ,  namely 


ip  (x)  = 

CO 

!  <|U  ~  I  a  •  (x)  .  (3.20) 

°  i=0  1 

This  recursion  scheme  is  consistent  with  the  formulation 
of  the  equations  for  cp  (x)  and  may  be  applied  directly. 
For  convenience,  define  U (x)  and  P(x)  by 


U  (x) 

=  k2(x) 

P  (x) 

(3.21) 

U*  (x) 

4U(x) 

(Note  that  the  discussion  here  is  valid  only  in  the 
allowed  region.  However,  by  making  the  mapping  (2.28), 
this  analysis  is  equally  valid  in  the  forbidden  region.) 
Then  observe  that 


a  (x)  = 

2ik  (x)  P(x) 

3(x)  = 

'  -  2ik(x)  P(x)  (3-22) 

Y  (x)  = 

1 

2ik (x) 

. 


, 
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By  the  explicit  use  of  the  preceding  recurrence  rela¬ 
tions,  one  finds 


ao(x)  =  2Tk7xT  p(x) 


Sc(x) 


Y  (x) 
'  o 


2ik (x) 


2ik (x) 


P  (x) 


(3.23) 


.  ,  ,  1  rU"  (X)  2  ,  \  n 

ao(x)  2ik (x)  {4U(x)  6P  (x)  }  * 


Noting  that 


P'  (x) 


U"  (x) 
4U(x) 


4P2 (x) 


(3.24) 


one  finds 


an  (x)  =  {  „ ,  h~,—  P  (x)  +  ^V7T7?-  -  3P2  (x)  }/C  (x) 


4ik (x) 


8U  (x) 


(x) 


=  ik  (x)P(x)/C  (x) 


Y i  (x) 


=  ik(x)/C(x) 


(3.25) 


a1  (x) 


{_  P4  (*)  +  P2  ( x )  P  '  (x)  +  U'"(x) 


2ik (x)  2ik (x) 


8U  (x) 


-  6P(x)P'  (x)  -  2  a-^(x)  [U1  (x) 


-  P  (x)P'  (x)]}/C(x) 


where  C (x)  is  defined  by 


‘ 
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C  (x)  =  2U  (x)  -  P2  (x)  . 


(3.26) 


Since  (x)  is  given  by 


{1  -  2a± (x) (x) } 


(3.27) 


cj)  (x)  may  be  approximated  by 


<J>  (x)  ~  aQ(x)  +  +  a2 


(3.28) 


Such  an  iterative  scheme  for  calculating  cf>  (x) 
can  be  carried  out  numerically  to  any  order  by  methods 
described  in  Chapter  IV  and,  where  it  shows  acceptable 
numerical  convergence,  leads  to  essentially  exact 
results  for  cf>  (x)  and  the  corresponding  solution  y  (x) 
to  the  Schrodinger  equation.  This  iterative  expansion 
for  (j)  (x)  appears  to  be  either  asymptotic  or  convergent 
and  is  very  much  more  accurate  than  the  JWKB  form  (1.3) 
over  essentially  the  entire  semiclassical  domain  where 
the  latter  provides  a  reasonable  approximation;  on  the 
other  hand,  when  the  JWKB  approximation  deteriorates, 
the  iterative  process  diverges. 

What  we  wish  to  emphasize  here,  however,  is  not 
the  detail  of  methods  for  evaluating  this  particular 
solution  4>  (x)  ,  but  rather  the  fact  that  such  a  particular 
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solution  exists  and  that  it  can  be  determined  to  extra¬ 
ordinary  accuracy  over  the  semiclassical  domain; 
furthermore,  as  noted  earlier  and  as  desired  originally, 
this  solution  is  smoothly  varying  and  does  not  present 
any  oscillatory  or  step-like  behaviour.  In  the  cases 
where  no  domain  of  semiclassical  behaviour  exists,  that 
is,  cases  for  which  the  JWKB  approximation  is  nowhere 
valid,  the  de  Broglie  wavelength  is  not  significantly 
smaller  than  the  characteristic  dimensions  of  the  po¬ 
tential  surface  and,  thus,  there  would  be  no  computa¬ 
tional  advantage  in  integrating  at  the  speed  of 
Hamilton's  equations  over  direct  numerical  integration 
of  the  second-order  Schrodinger  equation. 

Given  that  a  semiclassical  domain  exists  in  which 
an  accurate  "iterative"  estimate  of  this  particular 
solution  <J>  (x)  of  (2.13)  can  be  made,  we  must  then  con¬ 
sider  its  analytic  properties  so  that,  if  possible,  we 
can  obtain  its  approximate  analytic  or  numerical  con¬ 
tinuation  into  the  domains  where  the  iterative  scheme 
for  calculating  it  fails.  To  this  end,  we  now  consider 
properties  of  this  solution  (and  its  corresponding 
phase  correction  function  f  (x) )  in  the  turning  point 
region . 


' 


■ 
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D.  Behaviour  of  the  Desired  (x)  in  the  Turning  Point 


Region 


Since,  in  equations  (3.7),  the  roles  of  G(x)  and 
H(x)  are  interchanged  when  going  to  the  equation  for 
f  (x)  from  that  for  <j>  (x)  ,  it  is  tempting  to  consider  the 
use  of  the  iterative  solution  scheme  to  obtain  a  parti¬ 
cular  solution  for  f (x)  in  the  neighbourhood  of  a  turn- 

2 

ing  point,  where  k  (x)  =  0  : 


{1  -  f2  (x)  }  -  f'(x).  (3.29) 


In  Appendix  B,  we  discuss  the  particular  solutions  f  (x) 
which  can  thus  be  obtained.  However,  it  turns  out  that 
such  a  solution  for  f (x)  does  not  at  all  correspond 
with  the  particular  <J>  (x)  which  interests  us.  The 
asymptotic  behaviour  of  the  desired  <f>  (x)  in  the  neigh¬ 
bourhood  of  a  turning  point  may  be  deduced  as  follows: 

Recall  that  the  equation  which  determined  cj)  (x) 
is  (after  some  rearrangement) 


4> '  (x)  +  h(x)<j>(x)  =  g  (x)  { 1  -  (J)2(x)}. 


(3.30) 


Consider  the  substitution 


4)  (x)  =  Tanh  {U|X— } 


(3.31) 


The  differential  equation  then  obtained  for  u(x)  is 


- 


. 
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u'  (x) 


Sech‘ 


+  h(x)  Tanh 


{^} 


=  g (x)  Sech2  .  (3.32) 

In  the  vicinity  of  the  turning  point,  h(x)  vanishes 
and  one  must  solve 


u '  (x)  =  4g  (x) 


2k' (x) 
k(x) 


(3.33) 


This  is  an  exact  differential  and  one  finds  that 

=  in  fl M  (3.34) 

where  kQ  is  a  constant.  By  (3.31),  this  yields 

k (x)  -  k 

*(x)  =  kTxT~"k°'  (3'35) 

o 

and,  correspondingly, 

k 

f(x)  s  .  (3.36) 


Moreover,  this  result  for  the  phase  correction  function 
yields  fik  as  the  momentum  evaluated  at  the  turning 
point.  In  addition,  this  asymptotic  form  permits  <j)  (x) 
to  remain  within  the  unit  circle  and  provides  (j)  (x)  with 
a  limiting  value  of  -1  at  the  classical  turning  point. 
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However,  all  the  information  of  interest  to  us 
in  determining  the  particular  cf>  (x)  we  seek  is  in  the 
complex  constant  kQ  and  this  cannot  be  directly 
obtained  from  the  properties  of  $ (x)  in  the  asymptotic 
(semiclassical)  domain  but  requires  some  means  of 
propagating  the  exact  solution  across  the  intermediate 
region  to  the  turning  point.  By  showing  that  the  appro¬ 
priate  f (x)  has  the  form  (3.36)  near  the  turning  point, 
however,  we  have  definitely  excluded  the  two  types  of 
particular  solution  f  (x)  to  (3.29)  which  were  discussed 
in  Appendix  B. 

E.  Nature  of  the  Fundamental  and  General  Solutions 
(f)  (x)  and  f  (x) 

Additional  insight  into  the  structure  of  our 
problem  can  be  obtained  if  we  consider  the  collection 
of  the  various  fundamental  solutions  of  the  Riccati 
equations  (3.29)  and  (2.13)  (or  (2.31))  and  the  general 
solutions  which  can  be  constructed  from  these  funda¬ 
mental  ones. 

For  example,  if  we  examine  the  three  asymptotic 
forms  for  f (x)  in  the  turning  point  region  (namely  the 
two  obtained  in  Appendix  B  and  the  one  obtained  in  the 
preceding  section) ,  we  observe  that  they  are  distinct¬ 
ly  different  from  each  other.  The  reason  for  this  may 


■ 


« 
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be  traced  back  to  the  form  of  the  general  solution  to 
the  Riccati  equation  .  Given  a  function  ip(x)  de¬ 
termined  by 


ip  (x) 


g1 (x)  +  £g2 (x) 
g3 (x)  +  £g4 (x) 


(3.37) 


where  E,  is  a  constant,  and  where  g^ (x) ,  g2 (x) ,  g3  (x) 
and  g4 (x)  are  arbitrary  functions,  then  the  constant 
£  may  be  eliminated  by  obtaining  the  following  Riccati 
equation  for  ip  (x)  : 


i^'(x)  +  Q(x)^(x)  +  R(x)iJ;2(x)  =  P  (x)  (3.38) 

where 

Q  (x)  =  (g|(x)g4(x)  -  gx(x)g4(x)  +  g2(x)g3(x) 

-  g2  (x)  g3  (x)  }/D  (x) 

R (x)  =  (g3(x)g4(x)  -  g3 (x) g4 (x) }/D  (x)  (3.39) 

P(x)  =  (g|(x)g2(x)  -  g^ (x) g2 (x) }/D (x) 

D  (x)  =  g2(x)g3(x)  -  g1(x)g4(x) 


Although  four  arbitrary  functions  are  employed  in  de¬ 
fining  ip  (x)  ,  one  of  these  is  superfluous  since  (x) 
would  remain  unchanged  if  both  the  numerator  and 


denominator  of  (3.37)  were  to  be  divided  by,  for 
example,  g^ (x) .  Therefore,  given  four  solutions  to 
the  differential  equation  (3.38),  namely  (x)  , 

\l>2  (x)  f  ^3  (x)  an<3  ^4  (x)  (corresponding  with  constants 
^1'  ^2'  ^3  an<^  ^4  resPectivelY)  /  it  is  observed  that 
they  are  not  independent  since  their  "cross  ratio" 

(x)  -  (x)  }  {^2  (x)  -  ip4  (x)  }  (£x  -  K3)  (52  ”  £4) 

{ip1(x)  -  \p4  (x)  }{i|j2  (x)  -\p3(x)}~  (£ 1  -  C4)  (C2“  C3) 

(3.40) 

is  a  constant.  Thus,  three  solutions  to  the  Riccati 
equation  may  be  combined  to  provide  a  general  solution 
to  the  equation. 

We  have  found  two  asymptotic  forms  for  <p  (x)  that 

are  valid  in  the  JWKB  domain;  in  addition,  we  have  a 

technique  for  obtaining  these  to  arbitrary  accuracy. 

One  of  these  would  become  very  small  (going  as 
2 

k’(x)/4ik  (x) )  and  the  other  would  become  very  large 

2 

(going  as  4k  (x)/ik'  (x) ) .  Is  there  yet  a  third  funda¬ 
mental  solution  whose  magnitude  lies  between  the  modul 
of  the  other  two  solutions? 

Recall  that,  if  <£  (x)  satisfies  (2.13),  then 
1/4) *  (x)  also  satisfies  the  differential  equation. 
Hence,  <f>  (x)  may  possess  a  solution  on  the  unit  circle 
(since  4>  (x)  (x)  would  equal  unity).  Consider,  there¬ 

fore,  a  solution  of  the  form 


. 

■ 
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cf)(x)  =  -exp{-2i0  (x)  } 


(3.41) 


where  0 (x)  is  a  real-valued  function  which  satisfies 


0  '  (x)  =  k  (x)  +  Sin{2  0  (x)  }  . 


(3.42) 


Then,  the  wave  function  has  the  form 


y (x)  =  exp{ 


k (x) Cot [ 0  (x) ] dx}  . 


(3.43) 


Equation  (3.42)  may  be  solved  asymptotically  both 
far  away  from  and  in  the  vicinity  of  the  turning  point. 
When  k'(x)/2k(x)  is  much  smaller  than  k (x) ,  0 (x)  is 
given  asymptotically  by 


0  (x)  ~ 


k(x)  dx 


(3.44) 


x. 


Thus,  the  logarithmic  derivative  of  the  wave  function 
becomes 


x 


y1  (x) 
y  (x) 


55  k  (x)  Cot  { 


k (x)  dx} 


(3.45) 


*o 


The  latter  may  be  looked  upon  as  a  differential  equa¬ 
tion  and  has  the  immediate  solution 


x 


y  (x )  ~  Sin  { 


k(x)  dx} 


(3.46) 


x 
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Therefore,  this  third  fundamental  solution  to  the  equa¬ 
tion  for  <f>  (x)  provides  the  sinusoidal  form  that  charac¬ 
terizes  Calogero's  solution  and  the  real  JWKB  approxi- 
mants.  However,  in  obtaining  this  relationship,  we 
have  sacrificed  the  conceptual  link  between  quantum 
mechanics  and  classical  mechanics  since  it  is  no  longer 
possible  to  identify  momentum  with  a  continuous  function 
of  the  spatial  variable.  For  the  sake  of  completeness, 
consider  the  turning  point  region  where  equation  (3.42) 
becomes 


e'(x)  =  2k (x)~  sint2e(x)1  •  (3.47) 

This  equation  has  the  solution 

2  A  k  (x) 

0  (x)  =  Sin  {  ■  ^  -  "k'(x)  }  (3.48) 


where  kQ  is  a  constant,  and  has  the  corresponding  wave 
function 


x 


y  (x)  ~  exp {- 


k(x) 


o 


[k  -  k  (x) ]  dx}  . 
o 


(3.49) 


It  must  be  emphasized  that  this  is  a  formal  solution 
only  and  is  not  of  physical  interest.  Hence,  of  the 
three  fundamental  solutions  for  cp  (x)  ,  only  that  which 
may  be  obtained  by  the  recursive  scheme  described 
earlier  may  be  used  to  describe  the  motion  of  a  semi- 
classical  particle. 

*  *  * 
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The  first  chapter  had  presented  a  number  of  sug¬ 
gestive  arguments  for  considering  solutions  to  the  wave 
equation  in  complex  exponential  form.  However,  by  the 
end  of  the  second  chapter,  we  had  seemingly  deviated 
from  this  goal  and,  at  the  same  time,  introduced  the 
problem  of  solving  a  non-linear  equation  whose  behaviour 
was  more  obscure  than  that  of  the  Schrodinger  equation. 

In  this  chapter,  the  solution  was  cast  into  the  desired 
exponential  form  and  provided  momentum  eigenfunctions 
which  would  give  the  classical  momentum  in  the  limit  fi 
goes  to  zero.  We  identified  a  quantity,  namely  f (x) , 
as  a  correction  to  the  phase  and  amplitude  of  the  wave 
function.  The  properties  of  linear  fractional  trans¬ 
formations  and  of  the  differential  equations  for  f(x)  and 
c p  (x)  were  used  to  discuss  the  relationship  between  4>  (x) 
and  f (x) .  This  relationship  gave  rise  to  a  physically 
desirable  initial  condition  for  4)  (x)  and  provided  an 
upper  bound  to  its  magnitude  (Appendix  A) .  A  convergent 
iterative  scheme  for  evaluating  <j>  (x)  was  provided  such 
that  its  first-order  term  was  comparable  with  or  superior 
to  the  second-order  JWKB  approximant.  In  addition,  an 
asymptotic  form  for  cj>  (x)  valid  near  the  turning  point 
was  obtained.  Finally,  the  connection  between  this  so¬ 
lution  and  that  of  real-valued  wave  function  approxi- 
mants  was  shown.  Having  discussed  the  formal  properties 
of  this  technique,  it  is  now  necessary  to  consider  nume¬ 
rical  procedures  for  the  rapid  numerical  calculation  of 


the  wave  function. 
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CHAPTER  IV 

NUMERICAL  METHODS 


In  the  preceding  chapter,  a  method  for  solving 
a  particular  type  of  Riccati  equation  was  proposed: 
the  problem  was  reformulated  as  that  of  solving  the 
infinite  set  of  ancillary  equations  (3.16).  Since 
discussion  of  this  technique  has,  to  this  point,  been 
strictly  formal,  it  will  be  valuable  to  qualitatively 
consider  the  rationale  for  its  development. 

Consider  the  equation 

^Q(x)  =  aQ(x)  +  30(x)^q(x)  +  Y0(x)^0(x)  (4.1) 

where  aQ(x),  3q(x),  and  yQ  (x)  are  analytic,  small  com¬ 
pared  with  unity,  and  have  slowly  varying  derivatives 
for  some  domain  of  x.  Linearizing  this  equation  and 
neglecting  the  derivative  term  yields  aQ (x)  as  an 
approximate  solution.  Thus,  the  exact  solution  to 
this  equation  may  be  written  as 

ipQ(x)  =  aQ(x)  +  ^(x)  (4.2) 

where  ij^(x)  describes  the  correction  to  the  approximate 
solution.  Inserting  this  solution  into  equation  (4.1) 
results  in  a  Ricatti  equation  of  identical  form  to 
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(4.1)  where  the  coefficients  are  given  by  (3.18). 

Thus,  by  linearizing  the  solution  to  the  original 
Ricatti  equation,  we  obtain  a  new  Ricatti  equation 
whose  solution  describes  the  error  in  the  former. 

Since  the  task  of  solving  the  latter  equation  is 
formally  the  same  as  solving  the  equation  for  ^q(x), 
this  technique  may  be  iterated  (where  convergent) 
until  convergence  is  achieved. 

This  scheme  for  solving  a  non-linear  equation 
differs  markedly  from  the  JWKB  or  Langer  series.  The 
solution  to  all  orders  of  this  series  is  linearized; 
on  the  other  hand,  the  zeroth  order  term  in  the  Langer 
series  is  the  square  root  of  a  quantity  and,  ultimately, 
manifests  itself  in  the  Stokes  phenomenon.  In  addition, 
to  calculate  the  n-th  term  in  this  series,  we  need  only 
know  an_^  (x)  ,  3n_^(x),  an<^  Yn_i  (x)  /  the  JWKB  technique 
requires  that  all  preceding  terms  be  known  explicitly. 

As  a  consequence,  calculations  with  the  linearized 
iterative  scheme  are  less  laborious  and,  as  was  noted 
earlier,  more  accurate  than  JWKB  approximants  of  the 
same  order.  In  the  following  pages,  two  techniques 
for  pursuing  this  technique  are  described  and  tested. 
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A.  Rational  Function  Interpolation 


The  solution  to  ^q(x)  may  be  written 


^Q(x)  =  aQ  (x)  +  a-]_(x)  +  ...  +  aR_1(x)  +  ^(x),  (4.3) 


where  (x)  ,  i=0 , 1 ,2 , . . .n-1  are  known  and  where  ^n(x) 
satisfies  the  Ricatti  equation 


(4.4) 


Thus,  the  problem  of  solving  for  \pQ  (x)  is  reduced  to 
that  of  solving  for  ip  (x )  .  The  subscript  n  is  selected 
so  that  the  following  iterative  sequence  will  converge 
uniformly  in  some  domain  of  x: 


^n°} (x)  =  ao(x) 


(4.5) 


In  the  problem  of  calculating  cp  (x)  ,  n  appears  to  be 
unity  for  maximal  numerical  stability. 

In  order  to  utilize  equation  (4.5),  one  must 
select  a  finite  set  of  discrete  points  which  are  dis¬ 
tributed  in  some  regular  way  in  the  domain  of  interest. 
In  order  to  obtain  a  numerical  estimate  of  the  deriva¬ 
tive  on  a  set  of  discrete  points,  most  numerical 
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algorithms  will  use  finite  difference  estimates,  that 
is,  truncated  Taylor  series  expressions  of  low  order. 

In  our  problem,  however,  cj>  (x)  cannot  be  adequately 
described  by  a  polynomial  of  small  degree  since  the 
general  solution  to  a  Riccati  equation  necessarily 
has  the  form  of  a  ratio  of  two  functions.  To  overcome 
this  difficulty,  the  set  of  data  points  are  collocated 
with  a  "diagonal  rational  function  of  order  N , "  that  is 


rn(x) 


a,  +  a_x  +  ...  +  a  x 
1  2_ E 


P-1 


b.  +  b„x  +  ...  +  b  x 
12  q 


q-1 


N  =  p  +  q  -  1 


(4.6) 


q  <  p  <  q  +  1- 

(12) 

Rational  functions  are  widely  used  m  approximation 

theory  and  are  characterized  by  an  immunity  to  the  rapid 
fluctuation  of  derivative  that  often  accompanies  or¬ 
dinary  polynomial  interpolants .  These  features  make 
the  rational  function  ideal  for  this  iterative  procedure; 
moreover,  it  provides  a  simple  algebraic  form  for  (p  (x) 
instead  of  the  values  that  (f)  (x)  assumes  on  a  set  of 
points . 

For  convenience,  the  explicit  formulation  of  the 
iterative  procedure  for  the  allowed  region  is  provided 
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below.  For  its  application  to  the  forbidden  region, 
only  the  mapping  (2.28)  need  be  made. 

Recall  the  equation  for  the  function  $ (x) , 

namely 


(J)  (x) 


g(x) 

h  (x) 


U  -  <t>2  (X)  } 


4> 1  (x) 

h  (x) 


(4.7) 


where 


h  (x) 


g  (x) 


2ik (x) 

k  1  (x) 
2k  (x) 


(4.8) 


An  approximate  solution  to  (4.7),  which  shall  be  denoted 
<J) (x) ,  is  given  by 


4) 


(o) 


(x) 


g(x) 

h  (x) 


(4.9) 


Define  R(x),  which  is  the  residual  error  in  this  appro¬ 
ximation,  by 


R  (x) 


lh(x) 


[1  -  cj)(o)  (x)2] 


 i 


(O) 


(X) 


h  (x) 


<p(o)  (x)  }  (4.10) 


and  y (x) ,  which  is  the  correction  to  the  approximate 
solution  to  (4. 7), by 


y  (x)  =  <j>  (x)  -  (p  ^  (x) 


(4.11) 


f 
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The  resulting  equation 


y „(x)  =  {1  +  2  1  R(x) 

°  h  (x) 


(4.12) 


h  (x) 


i=0,l,2 


must  then  be  iterated  until  convergence  is  obtained. 

This  procedure  may  be  applied  to  both  bound 
state  and  scattering  problems;  the  only  significant 
difference  between  these  is  the  criterion  for  select¬ 
ing  the  discrete  points  at  which  cf>  (x)  is  explicitly 
evaluated.  For  bound  state  problems,  these  points  may 
be  equidistant;  this  is  also  the  case  for  part  of  the 
scattering  domain.  However,  in  the  region  where  the 
potential  "flattens"  (i.e.  asymptotically  tends  to 
zero  derivative  at  infinite  separation) ,  the  distri¬ 
bution  of  points  should  reflect  the  fact  that  the  more 
distant  portion  of  the  scattering  region  has  much  less 
effect  than  that  nearer  the  scatterer.  To  accommodate 
this  behaviour,  the  density  of  points  should  be  great¬ 
est  near  the  scatterer  and  least  in  the  region  where 
the  potential  flattens  out.  A  simple  way  to  assure 
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this  is  to  require  that  the  distance  between  successive 
pairs  of  points  decrease  proportionally.  For  example, 
one  might  let 


x.  .  ~  -  x .  , , 
i+2_ i+I 

x .  .  ,  -  x  . 
l+l  l 


=  r 


i=0 ,1,2 


/  •  •  • 


(4.13) 


where  T  is  a  constant  which  typically  lies  between  two 
and  five.  If  the  endpoints  xq  and  specified,  then 
x ^  must  be  selected  according  to 


x 


1 


(4.14) 


The  accuracy  of  this  procedure  depends  upon  two 
factors:  the  number  of  data  points  selected  and  the 

number  of  iterations  performed.  The  number  of  data 
points  employed  provides  an  upper  bound  to  the  accura¬ 
cy  of  this  technique.  The  reason  for  this  is  that  the 
number  of  data  points  corresponds  with  the  number  of 
rational  function  coefficients  and,  hence,  the  number 
of  variable  parameters  in  a  truncated  Taylor  series 
representation  of  $ (x) .  Typically,  ten  data  points 
will  provide  six  to  eight  significant  figures  in  an 
evaluation  of  <f>  (x)  .  From  numerical  investigations,  the 
rate  of  convergence  of  the  iterative  scheme  depends 
upon  the  Euclidean  norm  of  <j>  ^  (x)  for  the  domain  of  x 
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under  consideration.  Generally,  convergence  appears 
to  be  assured  if  the  domain  of  x,  namely  D,  is  selected 
such  that 

|  |  <j>  (x)  |  |  <  0.15  xeD  .  (4.15) 

For  a  fixed  number  of  data  points,  the  iteration  will 
converge  to  the  best  possible  approximation  to  (p  (x) 
for  the  number  of  variable  parameters  (that  is,  data 
points)  used.  In  order  to  estimate  when  convergence, 
in  this  case,  has  been  obtained,  one  need  only  iterate 
equations  (4.12)  until  the  Euclidean  norm  of  the  dif¬ 
ference  between  two  successive  approximations  to  y (x) 
is  less  than  some  tolerable  error.  It  is  important 
to  provide  an  upper  limit  to  the  number  of  iterations 
performed  in  a  computer  programme;  generally,  three  to 
ten  iterations  will  provide  estimates  of  <p  (x)  to  six 
digits  of  precision.  In  comparison  with  a  direct 
numerical  integration  of  Schrodinger ' s  equation  with 
six-figure  accuracy,  this  technique  provides  complex 
wave  functions  in  about  one-fourth  the  time.  If  <J>  (x) 
is  already  known,  the  calculation  of  the  wave  function 
is  an  order  of  magnitude  faster  than  a  direct  numerical 
integration  of  the  wave  equation . 

This  technique  was  the  first  used  in  order  to 
deduce  cj>  (x)  .  It  adds  to  the  repertory  of  numerical 
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analysis  a  new  technique  for  linearizing  and  iteratively 

solving  a  particular  form  of  the  Riccati  equation.  At 

the  same  time,  it  stimulated  the  belief  that  a  more 

simple  numerical  technique  with  a  greater  affinity  to 

analytical  methods  could  be  implemented.  Finally  and 

of  paramount  importance,  its  success  confirmed  that 
•  • 

Schrodinger ' s  equation  can  indeed  be  exactly  solved 
with  the  rapidity  of  Hamilton's  equation,  at  least  in 
the  semiclassical  region. 


B .  Linearized  Iterative  Sequence 

Consider  the  case  of  a  function  U(x)  which  is 
analytic,  single-valued  and  non-vanishing  in  some  domain 
D  (the  classically  allowed  or  the  classially  forbidden 
regions) .  Then,  this  function  possesses  a  local  Taylor 
series  about  every  point  in  D,  that  is  to  say  U(x) , 
where 

U  (x)  =  2m-<E--.-v-<*U 

tr 

=  k2 (x)  (4.16) 

may  be  associated  with  a  sequence  {LL,  i=0,l,2,...} 
where 


U±  =  U(x1) 


if  i=0 


dx' 


U(x) 


x=x. 


(4.17) 


if  i=l ,2,3,... 


X,£  D 
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By  making  similar  associations  of  derivative 
sequence  with  analytic  functions  a(x),  b(x),  c(x), 
the  so-called  "Leibnitz  rule  for  differentiation" 
of  a  product  may  be  written 

i 

c,  =  l  (h  a.  b.  .  i=0 ,1,2,...,  (4.18) 

j=0  3  3  3 

if  c  (x)  is  defined  by 

c(x)  =  a(x)b(x)  .  (4.19) 

Although  equations  (4.18)  and  (4.19)  are  equivalent, 
in  a  formal  sense,  they  are  very  different  to  compu¬ 
tational  devices  which  can  manipulate  only  numerical 
quantities,  but  never  algebraic  functions.  The 
Leibnitz  rule  can  be  extended  to  allow  for  many  other 
algebraic  operations. 

If  d(x)  is  the  first  derivative  of  c (x) ,  then 

di  =  ci+i  i=0,l,2,..  .  (4.20) 

2 

If  k  (x)  equals  U (x) ,  the  Leibnitz  rule  (4.18) 

if  i=0 

if  i=l  '  (4.21) 

i-1  . 

I  (L  k  .k ..}  if  i=2,3,4,... 
j=l  J  J  3 


may  be  written 


k. 

l 


±  JU~ 

v  o 


U. 


2k 


2K-{ui 
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The  indeterminacy  in  the  sign  is  due  to  the  selection 
of  the  branch  cut.  The  positive  case  corresponds 
with  a  left  going  wave,  the  negative  case  with  a 
right  going  wave. 

If  f(x)  is  defined  to  be  the  reciprocal  of  g  (x)  , 

then 


=  1/g 
1  ' 


if 


i=0 


i-1  . 

=  -f  y  (1)f.g.  . 

°  j£o  3  ^1-D 


(4.22) 


if  i=l ,2,3, 


Extensions  of  this  scheme  to  the  addition  of  two  func¬ 
tions  and  to  the  addition  or  multiplication  of  a  func¬ 
tion  by  a  constant  are  obvious. 

Upon  examination,  it  is  observed  that  the  linea¬ 
rized  iterative  scheme  described  in  the  previous 
chapter  may  be  employed  to  obtain  an  arbitrary  number 
of  terms  by  using  the  various  forms  of  Leibnitz  rule 
just  described  and  the  recursive  features  available  in 
modern  digital  computers.  Given  a  potential  where  the 
function  value  and  its  first  n  derivatives  are  known, 
a(x^)  and  its  first  n-1  derivatives  may  be  evaluated. 

In  calculating  each  term  in  the  a  series,  an  additional 
derivative  is  lost  so  that,  when  a  n  (xn )  is  evaluated, 
none  of  its  derivatives  are  known.  To  calculate  <}>(x^) 
to  arbitrary  accuracy,  it  is  only  necessary  to  take  a 
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"running  sum"  of  a^(x^)  function  values;  no  derivative 
information  need  be  stored.  This  technique  for  cal¬ 
culating  cj)(x^)  cannot  be  readily  extended  to  the  JWKB 
or  Langer  series  since,  in  order  to  calculate  the  n-th 
term,  the  preceding  n-1  terms  must  be  known  simulta¬ 
neously.  However,  as  each  term  in  the  a  series  relies 
only  on  the  immediately  preceding  one,  this  difficulty 
never  arises  and  it  is  not  necessary  to  store  or  mani¬ 
pulate  large  amounts  of  numerical  data.  It  must  be 
emphasized,  however,  that  this  technique  does  not 
provide  a  simple  algebraic  form  for  <p  (x)  (unlike  the 
method  of  rational  function  interpolants) .  Instead, 
it  can  provide  numerical  values  at  discrete  points  of 
almost  arbitrary  accuracy  at  specified  points.  If  an 
algebraic  form  for  cp  (x)  is  desirable,  the  numerical 
information  provided  by  this  method  can  be  collocated 
with  a  rational  function  or  other  interpolating  func¬ 
tion  . 

This  technique  has  proven  to  be  markedly  superior 
to  that  of  rational  function  interpolants.  It  is  sim¬ 
pler  to  use  and,  for  the  same  domain  of  x,  converges 
much  more  rapidly.  The  explicit  three-term  form  for 
cp  (x)  obtained  in  equations  (3.22)  to  (3.25)  will  give 
(J)  (x)  correct  to  seven  significant  digits  in  seventy  per 
cent  of  the  allowed  region;  four  or  five  digits  accuracy 
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is  available  in  eighty-five  per  cent  of  this  domain. 

It  is  vital  to  note  that,  as  the  three-term  series 
is  a  simple  algebraic  expression,  these  results  were 
obtained  with  only  a  modest  increase  in  computer  time 
compared  with  a  direct  integration  of  the  Hamilton 
characteristic  function.  The  ten-term  series  (ob¬ 
tained  from  the  application  of  the  Leibnitz  rule  to 
the  linearized  iterative  sequence)  would  converge  to 
ten-digit  accuracy  in  half  of  the  allowed  region; 
seven-digit  accuracy  has  been  obtained  for  seventy 
per  cent  of  the  region.  In  the  portion  of  the  allowed 
region  where  JWKB  approximants  deteriorated,  the  a 
series  oscillated  and  eventually  diverged. 

In  the  forbidden  region,  the  success  of  this 
technique  appeared  to  be  significantly  enhanced. 
Generally  speaking,  the  series  would  converge  when 
evaluated  at  the  same  distance  from  the  turning  ooint 
as  a  ooint  in  the  allowed  region  but  would  do  so  more 
rapidly.  In  fact,  convergence  of  ten  to  sixteen 
figure  accuracy  would  be  typical  and  would  be  obtained 
in  less  than  one-tenth  of  a  second.  Moreover,  the 
three-term  series  was  capable  of  providing  the  loga¬ 
rithmic  derivative  of  the  wave  function  with  ten  to 
twelve  digits  of  precision  when  evaluated  asymptotically 
far  from  the  turning  point.  It  is  essential  to  note 
that  no  known  techniques  for  numerical  integration 
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reliably  predict  the  logarithmic  derivative  to  more 
than  six  or  eight  digits  of  accuracy.  Thus,  the 
linearized  iterative  technique  (or  its  explicit  three- 
term  form)  has  realized  the  objective  of  numerically 
calculating  the  wave  function  with  the  ease  and 
rapidity  one  would  associate  only  with  the  solution 
to  Hamilton's  equations. 

By  the  iterative  process  described  in  this 
chapter,  we  may  obtain  essentially  exact  values  for 
the  desired  function  cp  (x)  and,  hence,  for  the  corres¬ 
ponding  exact  "pure  complex  exponential"  solution  to 
the  Schrodinger  equation  in  the  region  where  a  semi- 
classical  (JWKB)  approximation  would  be  valid.  We 
have  given  some  brief  exemplary  data  on  numerical  con¬ 
vergence  of  the  iterative  process  in  the  two  cases. 

Again,  we  wish  to  emphasize  that  the  objective 
of  this  thesis  is  not  to  present  a  complete  set  of 
numerical  techniques  for  calculating  the  solution 
y(x)  to  Schrodinger ' s  equation  but  to  show  that,  in 
principle,  the  exact  wave  function  may  be  obtained 
everywhere  at  a  speed  controlled  by  Hamilton's  equations 
and  that  this  is  practicable  in  the  region  of  rapid 
oscillation  in  the  Schrodinger  wave  function  if  we 
preserve  the  conceptual  link  between  quantum  and 
classical  mechanics  in  the  method  employed.  This  we 
have  done  by  demonstrating  the  existence  of  a  smooth 
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wave  mechanical  correction  function  <j>  (x)  and  calcula¬ 
ting  it  "iteratively"  in  the  semiclassical  regions. 
However,  for  the  sake  of  completeness,  we  shall 
briefly  discuss  the  problem  of  rapid  numerical  con¬ 
tinuation  through  the  turning  point  region  and  neigh¬ 
bouring  domains  where  the  iterative  scheme  for  <p (x) 
fails . 

We  may  begin  by  recalling  that,  in  the  classi¬ 
cally  forbidden  regions,  the  solution  y(x)  =  ip  (x)  , 
say,  generated  by  the  cf>  (x )  iterative  scheme  is  just 
the  regular  solution  which  is  normally  of  interest. 

In  the  classically  allowed  region,  the  iterative  solu¬ 
tion  (J)  (x)  to  (2.13)  and  its  complex  conjugate  <f>*  (x) 
generate  a  pair  of  complex  exponential  solutions  y(x), 
y* (x) ,  and  some  (conventionally  real)  linear  combina¬ 
tion  of  these , 

^ (x)  =  ay (x)  +  a*y*  (x)  ,  (4.23) 

which  is  the  analytic  continuation  of  the  solution 
originated  in  a  forbidden  region.  To  determine  the 
coefficient  a  requires  some  process  for  propagating 
ip  (x)  across  the  region  between  the  two  domains  where 
iterative  solution  of  (2.31)  or  (2.13)  can  be  done. 

This  region  has  essentially  two  distinct  parts 
which  must  be  treated  differently:  (a)  the  turning 


’ 


' 

a 


' 


61 


point  region  proper,  where  |k  (x) |  <  |k*  (x) | ,  and 
(x)  has  a  form  described  approximately  by  (3.35); 

(b)  an  intermediate  region  where  |k  (x)  |  ~  | k  *  (x)  | 
and  the  iterative  device  converges  poorly.  In  the 
turning  point  region  proper,  a  variety  of  techniques 
may  be  considered,  including,  for  example,  methods 
based  on  various  power-series  solutions.  These  will 
suffice  to  join  a  solution  ip  (x)  coming  from  the  conver¬ 
gence  region  of  the  iterative  device  (2.31)  or  (2.34) 

in  the  non-classical  region  and  carry  it  across  the 
2  ' 

zero  of  k  (x) ;  however,  they  deteriorate  badly  as  soon 
as  the  oscillations  of  ip  (x)  commence  in  the  classically 
allowed  region  and,  using  them,  we  can  at  best  propa¬ 
gate  ip  (x)  to  its  first  node.  Between  this  point  and 
the  semiclassical  domain  where  the  iterative  solution 
for  <p  ( x )  is  known  lies  an  "intermediate"  region.  In 
this  region,  several  oscillations  of  ^ (x)  maY  occur  so 
that  direct  numerical  integration  of  the  Schrodinger 
equation  would  be  too  time-consuming.  However,  in  this 
intermediate  region,  numerical  integration  of  (2.13) 
or  some  smooth  related  differential  equation  will 
become  stable  as  the  asymptotic  iterative  method  becomes 


unsuitable . 


, 
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CHAPTER  V 

CONCLUSION 

In  the  preceding  chapter,  we  investigated  a  set 
of  techniques  for  isolating  a  particular  solution  to 
the  equation  for  $ (x) ,  one  which  asymptotically  pro¬ 
vided  the  JWKB  approximants .  We  know  that  any  solution 
to  the  equation  for  <f>  (x)  furnishes  a  pair  of  solutions 
to  the  one-dimensional  Schrodinger  equation  and  yet 
there  is  an  especially  attractive  quality  surrounding 
this  particular  choice  of  <j>  (x)  . 

In  the  appendix,  we  observe  that  <j>  (x)  is  invari¬ 
ably  either  inside  the  unit  circle,  outside  the  unit 
circle,  or  on  the  unit  circle.  In  terms  of  the  linear 
fractional  transformation  which  relates  <p  (x)  to  f  (x)  , 
the  phase  correction  function  is  either  in  the  right- 
half  complex  plane,  in  the  left-half  complex  plane,  or 
on  the  imaginary  axis,  respectively.  Therefore,  the 
theorem  and  corollaries  of  Appendix  A  may  be  looked 
upon  as  a  form  of  conservation  principle,  that  is  to 
say,  the  preservation  of  the  direction  of  a  quantum 
mechanical  particle's  momentum.  If,  initially,  a  par¬ 
ticle  may  be  described  as  being  predominantly  left 
going  (  |  cj)  (x  )  |  <  1)  ,  it  is  always  so.  If,  initially, 
a  particle  may  be  described  as  being  predominantly 
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right  going  (|<j>(xo)|  >  1),  it  is  always  so.  Finally, 

if  a  particle  may  be  described  as  possessing  equi- 
probable  left  and  right  going  states  (  |  (j)  (x  )  |  =  1)  , 
it  will  always  remain  so.  Of  all  the  left  going 
solutions,  the  techniques  of  the  preceding  chapter 
provides  one  which  is  uniquely  that  of  a  pure  left 
going  wave,  an  observation  which  was  verified  by  com¬ 
parison  with  a  direct  numerical  integration  of 
Schrodinger 1 s  equation.  What  then  may  one  conclude 
about  other  solutions  to  <j)  (x)  which  are  contained 
inside  the  unit  circle? 

In  the  third  chapter,  we  saw  that  a  general 
solution  to  a  Riccati  equation  could  be  obtained 
from  three  particular  solutions.  In  the  case  of  the 
equation  for  cf>  (x)  ,  one  might  select,  as  particular 
solutions,  the  one  obtained  from  the  linearized  itera¬ 
tive  sequence,  the  complex  conjugate  of  its  reciprocal, 
and  one  which  is  of  unit  modulus.  It  is  important  to 
note  that  the  latter  solution  is  characterized  by  a 
double  wave-frequency  clockwise  rotation  on  the  unit 
circle  (see  equations  (3.49)  and  (3.50)).  On  the  other 
hand,  the  solution  to  (j>(x)  which  we  have  sought  des¬ 
cribes  a  contour  which,  commencing  with  a  point  on  the 
negative  real-axis  (corresponding  with  the  zero-force 
point),  proceeds  slowly  upward,  parallel  to  the 
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imaginary  axis.  Then,  it  follows  a  sweeping  retrograde 
motion  that  goes  through  180°  becoming  asymptotically 
parallel  to  the  imaginary  axis  as  the  classical  turn¬ 
ing  point  looms  nearer.  All  other  solutions  to  the 
equation  for  c p  (x)  which  are  bounded  by  the  unit  circle 
contain  a  "mixture"  of  these  three  fundamental  solutions 
(owing  to  the  linear  fractional  transformation  involved) . 
This  has  two  deleterious  effects.  It  contains  some  of 
the  character  of  the  unit  modulus  solution,  that  is, 
the  double  frequency  oscillations  which  have  plagued 
many  numerical  investigations  of  the  solution  to 
Schrodinger 1 s  equation.  At  the  same  time,  this  solu¬ 
tion  does  not  provide  a  pure  left  going  wave  but  con¬ 
tains  a  limited  amount  of  right  going  character.  By 
analogy  to  classical  wave  motion,  this  corresponds  to 
a  reflected  wave  which  manifests  itself  as  a  form  of 
interference  and  is  typified  by  a  double  local  wave- 
frequency.  This,  then,  has  the  effect  of  introducing 
spurious  fluctuations  in  the  real  and  imaginary  parts 
of  the  logarithmic  derivative  of  the  wave  function. 

Thus,  the  deceptive  undulatory  appearance  of  numerical 
estimates  of  the  quantum  mechanical  momentum  are  noth¬ 
ing  more  than  the  result  of  not  having  selected  the 
correct  initial  condition  for  a  pure  travelling  wave. 

The  linearized  iterative  sequence  provides  a  pure  left 
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or  right  going  wave  and,  hence,  a  quantum  mechanical 
momentum  which  departs  from  its  classical  counterpart 
in  a  smoothly-varying,  non-oscillatory  manner.  More¬ 
over,  it  does  so  in  a  manner  which  permits  the  calcu¬ 
lation  of  the  wave  function  with  the  ease  and  rapidity 
of  solving  the  Hamilton  equations  of  motion. 

The  strength  of  this  technique  is  that  it  embodies 
a  methodological  link  between  the  wave  and  particle 
descriptions  of  matter.  By  reformulating  the  problem 
of  solving  the  one-dimensional  Schrodinger  equation 
as  a  correction  to  an  asymptotic  form,  one  preserves 
the  intimate  connection  between  the  local  and  wave 
mechanical  momenta.  By  confining  onself  to  a  parti¬ 
cular  solution  of  the  resultant  Riccati  equation,  the 
wave  function  necessarily  becomes  the  exponential  of 
a  modified  form  of  the  Hamilton  characteristic  function. 
Simultaneously,  the  wave  function  becomes  numerically , 
no  more  elusive  than  the  solution  to  the  Hamilton- Jacobi 
equations . 
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APPENDIX  A 


BOUNDING  OF  <J>  (x) 

Theorem: 

If  1 4>  ( xq )  |  <  1  for  some  xq  in  the  classically 
allowed  region  D,  then  <j>  (x)  is  bounded  by  the  unit 
circle  for  all  x  in  D. 

Proof : 

The  function  <f>  (x)  is  defined  by 


<f>  (x) 


U-  *2(x)} 

4ikz (x) 


4> 1  (x) 

2ik (x) 


(A .  1 ) 


To  <j)  (x)  associate  a  real-valued  non-negative  modulus 
r(x)  and  a  real-valued  phase  0 (x)  such  that 


<J>  (x)  =  r(x)  exp{i0(x)} 


(A. 2) 


Coupled  equations  in  r (x)  and  0 (x)  arise,  namely 


r*  (x) 


k*  (x)  r, 

2k(x)  1 


r2  (x) }  Cos{ 0 (x) } 


0' (x)  =  -2k (x) 


k'  (x) 
2k  (x) 


Sin [0 (x) ] 


1+r2 (x) 

r  (x) 


(A. 3) 


Given  the  initial  condition 


r  =  r  (x  )  <  1 
o  o 


(A. 4) 
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a  general  solution  to  the  first  equation  in  (A. 3)  is 


r(x)  =  exP  m  — 
exp  u(x)  +  c 


(A. 5) 


where  u(x)  is  given  by 


u(x)  =  2 


k'  (x) 
2k  (x) 


Cos { 0  (x) }  dx 


x 


(A. 6) 


and  where  c  is  a  real-valued  constant;  the  initial  con¬ 
dition  confines  c  to  the  value 


c  = 


1  -  r 

c 

1  +  r 


(A. 7) 


Since  rQ  <  1 ,  c  must  be  positive.  For  an  arbitrary 
value  of  xeD,  (A. 5)  may  be  rewritten  as 


1  -  r(x)  =  c{l+r(x)}  exp{-u(x)  } 


(A. 8) 


But,  the  right-hand  side  of  the  latter  is  positive 
definite  and 

1  -  r (x)  >  0  (A. 9) 

completing  the  proof. 

*  *  * 


Corollary : 

If  |  cf)  (x  )  |  >  1  for  some  xq£:D,  then  |  <{>  (x)  |  >  1 


for  all  x  e  D. 
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Proof ; 

A 

Define  a  function  4> (x)  by 


$  (x)  =  1/4) *  (x)  .  (A. 10) 

A 

By  direct  substitution,  it  is  observed  that  41  (x)  satis¬ 
fies  (A.l)  and  that  |4>(xq)  |  <  1.  By  the  preceding 
theorem,  |  $  (x)  |  <  1  for  all  xeD.  Thus,  1 4>  (x)  |  >1  for 
all  xeD. 

*  *  * 


Corollary : 

If  [<j)(x  )  |  =  1  for  some  xqsD,  then  |4>(x)  |  =  1 
for  all  xeD. 

Proof : 

Write  4>  (x)  in  polar  form  (A.  2);  then,  r(x)  is 
given  by  (A. 5)  where,  by  (A. 7),  c  is  identically  zero. 
Therefore,  from  (A. 5), 

r(x)  =  1  (A. 11) 


for  all  xsD. 
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APPENDIX  B 


TWO  ASYMPTOTIC  FORMS  FOR  f (x) 


In  calculating  f (x)  at  the  turning  point,  one 
would  prefer  to  employ  the  iterative  scheme  on  the 
differential  equation  for  the  phase  correction  func¬ 
tion.  Unfortunately,  the  derivative  of  a(x)  for  this 
case  is  not  slowly  varying.  So,  for  the  solution  to 


(B.l) 


it  is  not  sufficient  to  say  that  an  asymptotic  solution 
to  f (x)  valid  in  a  neighbourhood  about  the  turning 


2 

point  is  ik  (x)/k' (x) ,  since  the  derivative  of  this 


quantity  is  not  slowly  varying.  A  more  flexible  asymp- 


2 

totic  form  is  cik  (x)/k'  (x)  where  c  is  a  constant  that 


is  to  be  determined.  Detailed  calculations  show  c  to 
be  1/4  and  that 


(B .  2 ) 


The  corresponding  wave  function  is 


y (x)  =  exp  {-  j 


(B .  3 ) 


and,  in  the  case  of  a  linear  turning  point, 
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k2 (x)  =  x  (B. 4) 

this  becomes 

y(x)  s:  exp[-x3/6]  .  (B.5) 

This  solution  is  related  to  the  Airy  functions  by 

B. (-x)  r(i)  31/3 

y  (x)  =  {A.  (-x)  +  - }  - 5 -  .  (B.6) 

1  /3  2 


Recall  that,  in  the  case  of  the  differential  equation 

for  <p  (x)  ,  the  function  l/(f>*  (x)  is  also  a  solution. 

2 

This  suggests  that  some  constant  times  k' (x)/ik  (x) 
might  be  an  asymptotic  solution  for  f  (x) .  Similar 
calculations  show  that  f (x)  may  be  given  approximately 
by 

f  (x )  «  -  2l^'  (X-^-  (B .  7 ) 

kZ(x) 


and  that  the  wave  function  becomes 


x 

y (x)  =  exp  { 


d  n  2  ,  « .  , 

(x  )] 

,2  ,  i  v 

k  (x  ) 


dx'  }  . 


(B .  8 ) 


In  the  case  of  a  linear  turning  point,  the  latter 
becomes 


y  (x)  ~  x 


(B .  9 ) 


which  is  related  to- the  Airy  functions  by 


. 
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y (x)  S  (Ai(-x)  - 


r(j)  31/3 


(B.10) 


/3 


Unfortunately,  these  asymptotic  forms  do  not 
conform  to  the  type  of  solution  one  would  require  if 
it  were  to  be,  at  the  same  time,  a  momentum  operator 
eigenfunction.  The  momentum  of  a  semi-classical 
particle  "tunneling"  through  the  turning  point  into 
the  forbidden  region  should  have  a  finite  but  non¬ 
vanishing  "momentum."  Therefore,  the  solutions  ob¬ 
tained  here  for  f (x)  must  be  discarded  in  favour  of 
one  which  behaves  like  some  constant  terms  l/k(x)  in 
the  vicinity  of  the  turning  point. 
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APPENDIX  C 

SOME  ILLUSTRATIVE  NUMERICAL  RESULTS 

Having  made  a  number  of  comments  regarding  the 
properties  of  the  linearized  iterative  sequence,  we 
give  here  some  illustrative  examples  of  numerical 
results  which  demonstrate  the  ease  and  accuracy  in¬ 
herent  to  this  technique.  As  a  model,  the  twenty- 
fourth  level  of  a  simple  harmonic  oscillator  potential 
is  tested  below.  There  is,  however,  no  reason  to 
suspect  that  these  results  are  peculiar  to  the  harmonic 

oscillator;  in  fact,  this  problem  is  a  "worst  possible 

2 

case"  since  [k'(x)/k  (x) ]  does  not  go  to  zero  as 
strongly  for  this  potential  as  for  most  physically 
realistic  potentials. 

*  *  * 

Classically  Forbidden  Region 

In  the  case  of  the  forbidden  region,  we  numeri¬ 
cally  integrate  the  Schrodinger  equation  with  a  fourth 
order  Runge  Kutta  from  15.0  to  the  turning  point  (7.0). 
At  selected  points  (denoted  X) ,  we  output  the  logari¬ 
thmic  derivative  (YPRIME/Y)  from  the  exact  numerical 
integration  (EXACT) ,  from  a  second  order  JWKB  calcu¬ 
lation  (WKB)  and  from  a  three-term  a  series  for 


' 
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(SERIES).  Similarly,  the  exact  and  series  estimates 
of  f  (x)  (denoted  F)  and  cj)  (x)  (denoted  PHI)  are  given. 
For  example,  at  14.9  (which  is  "far"  from  the  turning 
point) ,  we  have 

X 

YPRIME/Y (EXACT) 

YPRIME/Y (WKB) 

YPRIME/Y (SERIES) 

F  (EXACT) 

F (SERIES) 

PHI (EXACT)  -1.625783097676806D-03 

PHI (SERIES)  -1. 6 2516 43 3520 770 8D-0 3 

Consider,  also,  the  following  estimate  for  <J>  (x)  in 
terms  of  a  nine-term  a  series.  Note  the  extraordinary 
rate  of  convergence  of  the  series. 


1. 489999999999999D  01 
-1.319616513368636D  01 
-1.319638767016770D  01 
-1.3196148803069 8 7D  01 
1. 00325686114515 4D  00 
1. 00325561958725 7D  00 


PHI  -1.62516036 87 D- 03 

SERIES 


-1. 636889895 8D-03 
1.190470491 2D-05 
-1.7916186870D-07 
4 . 10  5  82  76  33  8D-0  9 
-1.2643232256D-10 

4 . 888222  845 5D-1 2 
-2. 274140561 6D-13 
1 . 2  36  55  816  91D-14 
-7.6938699336D-16 


1 
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By  direct  comparison,  it  may  be  noted  that  the  three- 
term  a  series  estimate  of  (J)  (x)  is  more  accurate  than 
the  exact  numerical  integration!  Indeed,  far  from 
the  turning  point,  the  three- term  series  will  typically 
provide  ten  digit  accuracy. 

Now,  let  us  look  at  some  results  obtained  signi¬ 
ficantly  closer  to  the  turning  point  but  still  in  the 
asymptotic  domain,  namely 


X 


1. 09999999999999 9D  01 


YPRIME/Y (EXACT)  -8 . 560416092333452D  00 


YPRIME/Y (WKB) 


-8. 56167026 3127 45 3D  00 


YPRIME/Y (SERIES)  -8 . 560420804163676D  00 


F  (EXACT) 


1.008854711444572D  00 


F (SERIES) 


1.008855266739089D  00 


PHI (EXACT) 


-4. 407840643788936 D- 03 


PHI (SERIES) 


-4. 408115848716355D-03 


and 


. 
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PHI  -4.4078401422D-03 

SERIES 


-4 . 5012  5  84  451D-0  3 


9 . 7542560002D-05 


■4.4011171042D-06 


3. 0178963824D-07 


-2 . 7741564162D-0  8 


3. 196122616 3D-09 
-4 . 4252879354D-10 
7 . 1548228631D-11 
-1.3228271295D-11 


At  this  point,  the  rate  of  convergence  of  the  a  series 
is  no  longer  so  rapid  and  the  estimate  of  (p  obtained 
from  the  three-term  a  series  is  less  precise  than  that 
obtained  from  the  numerical  integration.  However,  it 
should  be  noted  that  the  three-term  a  series  evaluation 
of  (f>  (x)  is  correct  to  the  seventh  significant  figure 
and,  thus,  is  sufficiently  accurate  for  most  physical 
applications . 

Consider,  finally,  the  application  of  this  tech¬ 
nique  to  a  point  in  the  near  vicinity  of  the  turning 
point.  Here,  the  a  series  is  asymptotic  and  begins  to 
diverge  after  the  sixth  term,  namely 


■ 
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X 

YPRIME/Y (EXACT) 
YPRIME/Y (WKB) 
YPRIME/Y (SERIES) 
F (EXACT) 

F (SERIES) 

PHI (EXACT) 

PHI (SERIES) 

and 

PHI 

SERIES 


7. 99999999999999 3D  00 
-4.108657178764820D  00 
-4. 13965001287407 ID  00 
-4.115851259271729D  00 
1.060850721908783D  00 
1.062708225508418D  00 

-2. 9 526 9 9 15 7 77416 6D-0 2 
-3 .0400918914725 82D-02 

-3.1253411577D-02 


-3. 442651863 3D-02 
6 . 5807585204D-03 
-2 . 5600978225D-03 
1. 4974732945D-03 
-1. 1671202290D-03 

1. 1356912060D-03 
-1 . 32423844 47D-03 
1.7987116419D-03 
-2. 788071110 8D-03 


It  was  observed  that,  for  the  case  of  the  twenty-fourth 
harmonic  oscillator  level,  the  three-term  a  series  was 
sufficiently  accurate  in  estimating  <j)  (x)  for  x  greater 
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than  9.0  to  satisfy  any  reasonable  computational 
requirements  in  calculating  the  wave  function  in 
the  forbidden  region,  at  the  same  time,  provided  the 
solution  with  the  speed  of  solving  Hamilton's  equations. 

*  *  * 

Classically  Allowed  Region 

The  a  series  is  rapidly  convergent  for  seventy 
per  cent  of  the  allowed  domain.  Fifteen-term  a  series 
are  given  below  for  a  set  of  points.  (Note  that  the 
numerical  values  for  <p  and  a.  have  two  components  cor¬ 
responding  with  real  and  imaginary  parts.) 


X 

o 

• 

o 

PHI 

-5. 200 7 6 149 8 9 336 4 OD-O 5 

o 

• 

o 

SERIES 

1: 

o 

• 

o 

o 

• 

o 

2: 

-5. 20616409 8292376D-05 

0.0 

3: 

0.0 

o 

• 

o 

4: 

5. 420 89 9 4 779 8200 9D-0 8 

o 

• 

o 

5: 

8. 4 577017470 13455D-13 

o 

• 

o 

6: 

-1. 85143771245299 8D-10 

o 

• 

o 

7: 

-8. 79 8747020 599 331D-15 

o 

• 

o 

8: 

1.321319924872674D-12 

o 

• 

o 

9: 

1.271007610222153D-16 

0.0 

10: 

■  -1.612967304265969D-14 

o 

• 

o 

11: 

-2.609962914660366D-18 

o 

• 

o 

12: 

3. 00 45 87 879 310 2 3 8D-16 

o 

• 

o 

13: 

7. 3450 87590 179455D-20 

o 

• 

o 

14  : 

-7.931679272439819D-18 

o 

t 

o 

15: 

-2 .72 9 65 56 03 1716 9 0D- 21 

o 

• 

o 

■ 


X 


5.000000000000000D-01 
-5. 334 7540 54 3 86 6 51D-0 5 
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It  has  been  observed,  heuristically ,  that  the 
allowed  region  may  be  considered  to  have  three  dis¬ 
tinct  regions.  First,  there  is  the  region  where  the 
a  series  effectively  converges.  This  is  approximately 
seventy  per  cent  of  the  allowed  region  and  the  three- 
term  series  for  cj>  (x)  will  give  seven-figure  wave 
functions  in  most  of  this  domain.  Then,  there  is  a 
region  (about  twenty  per  cent  of  the  allowed  domain) 
where  the  a  series  is  asymptotic.  Finally,  the  re¬ 
maining  ten  per  cent  of  the  domain  (namely  that 
localized  about  the  turning  points)  is  completely 
divergent.  However,  it  must  be  emphasized  that  it 
is  in  the  region  where  the  linearized  iterative 
sequence  is  successful  that  numerical  integration 
techniques  are  most  impaired  and  that  it  is  here  that 
we  can  readily  solve  the  one-dimensional  Schrodinger 
equation  with  the  speed  of  Hamilton's  equation.  More¬ 
over,  we  have  already  established  stable  numerical 
integration  techniques  which  allow  us  to  extend  the 
solution  to  the  turning  point. 
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